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Abstract— This paper deals with the buckling of the stiffened plate under uni-axial compression. The 
direct integration of the governing differential equation is performed and the exact solution to the 
problem is obtained. As examples, a square plate with single stiffener, and a stiffened three-span, 
continuous plate are investigated, with special attention given to the influence of stiffener misplace- 
ment on the buckling load and mode shape of the plate. It is found that a small misplacement of the 
stiffeners from the nominal configuration may change the buckling mode from a global one to a 
highly localized one. 


1. INTRODUCTION 

Traditionally, the stability of the stiffened plate has been studied following three different 
trends. One approach consists of replacing the stiffened plate by an ‘equivalent’ orthotropic 
plate after the stiffeners are smeared out, in an energetic sense, over the entire surface of 
the plate [1, 2]. This approach appears reasonable for plates with many closely spaced 
stiffeners but is doubtful for plates with fewer stiffeners. The second approach is based on 
energy consideration and treats the contributions of the plate and the stiffener separately; 
the Rayleigh-Ritz method has been utilized widely to estimate the buckling load of the 
stiffened plate structure [3, 4]. This method may predict well the global buckling but fail to 
detect the localization of buckling mode due to the small changes in the location of the 
stiffeners. The third approach is the analytical method for equally spaced stiffeners by the 
analytical finite difference calculus [5, 6] which, though powerful for studying plates with 
periodically spaced stiffeners or supports, is totally inapplicable if the periodicity is 
disturbed as is usually the case when misplacement in the location of the stiffener or 
support is present through imprecision of construction. Despite their usefulness and 
simplicity, the above-mentioned methods can only be employed to investigate the global 
buckling of the structure and appear incapable of revealing the localization phenomena 
when the structure is sparsely or irregularly stiffened and the buckling mode is likely to be 
localized. 

In this paper, we investigate the effect of small structural irregularity, due to the 
misplacement of stiffeners or interior supports, on both the buckling load and the buckling 
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mode of the rib-stiffened plate. Since the buckling mode shape is of main interest here, the 
interaction between the plate and stiffeners should be properly taken into account; and 
none of the above-mentioned methods seems successful in this context. Here, the 
integration of the general governing differential equation is attempted for the stiffened 
elastic plate. By considering the rib-stiffened plate as a physically continuous plate with as 
many spans as the number of ribs, the stiffeners are accounted for through the conditions 
of continuity. Two cases commonly encountered in practice are considered; one with 
simple support under the ribs and one without. It is found that in the presence of small 
misplacement of stiffeners or interior supports, the buckling mode shapes experience 
dramatic changes to become strongly localized. Localization phenomenon was first un- 
covered by the Nobel Laureate P. W. Anderson [7] in physics. Its occurrence in structures 
has recently attracted much attention. Among others, Pierre and Plaut [8] considered the 
two-span column case with deterministic disorder. A more general case, the multi-span 
column, was recently treated by Nayfeh and Hawwa [9] using the transfer matrix method. 
Ariaratnam and Xie [10] investigated the localization in the buckling of a system of rigid 
bars connected with springs in the stochastic setting. Tvergaard and Needleman [11] 
discussed the development of localized patterns in the elastic-plastic and thermal buckling 
problems. Cai and Lin [12] studied the localization of wave propagation in randomly 
disordered periodic structures. The deterministic buckling localization in cylindrical shells 
was investigated by El Naschie [13-16]. In this study, we deal with the localization 
phenomenon in the buckling of stiffened plates. As a numerical example, a two-span plate 
with a single rib is discussed using different parameters for the stiffener. Furthermore, a 
stiffened three-span plate is also investigated, and the optimal configuration of stiffener 
placement, which yields the highest buckling strength, is discussed along with the attendant 
localization sensitivity to deterministic misplacements. 


2. STABILITY FORMULATION 


We consider a rib-stiffened rectangular plate subjected in its mid-plane to uniform 
compression P in the x direction (Fig. 1). The differential equation of the deflection 
surface of the plate under consideration is 


d {^ + 2 -^_ + £2L) + - o, (l) 

\ 3x 4 3x 2 3_y 2 3y 4 / dx 2 

where w is the transverse displacement, downward positive; D is the flexural rigidity of the 





Fig. 1. Simply supported rectangular plate stiffened with N ribs. 
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plate. The solution of equation (1) can be represented in the following form 

w(*, v) = X(x)sin|-^-j. 

Substitution of equation (2) into equation (1) results in 

*X + (L- 2^ + ^X = 0. 

dx 4 \D b 2 / dx 2 b 4 


The corresponding characteristic equation reads 


or 



( 2 ) 

(3) 

(4) 

(5) 


Even for the unstiffened plate, the buckling load P cr is always equal to or larger than 
A-^D/b 2 [17]. Thus, for rib-stiffened, or intermediately supported, plates, we have roots 
s,(i = 1, 2, 3, 4) as following 

Sj = ifa, s 2 = -ifa, s 3 = ifa, *4 = ~ifa, (6) 


where 





Solution of equation (1) can be written as 


w(r) = [,4 cos (ft*) + B sin (fax) + Ceos (fax) + D sin (fax)] sin 



(7) 

( 8 ) 


where A, B, C and D are unknown constants, which are to be determined by use of 
continuity and boundary conditions. For the arbitrary, ;'th span, the solution can be written 
as 

Wj(Xj) = [Aj cos ( faxj) + Bj sin (faxj) + C, cos (fax,) + D f sin (ft*,)] sin 

(9) 



where a ; is the length of the jth span and j ranges from 1 to TV for an TV-span plate. We 
consider the plate simply supported along its periphery. Then the boundary conditions are 


w 


llxj* 


=o — 0, 


- -»(-^T + *— p) 
\ 3>Vi 3 y 2 I 

r W| = - D (^ + v ^\ 

XN “ N \ 3 x 2 n 3/ / 


= 0, 


*>=o 


= 0, 


x N~ a N 


Wjvl 


x N~ a N 


= 0, 


( 10 ) 
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where Af* 1 ’ and M[ N) are the bending moments in the first and last spans of the continuous 


plate; v is the Poisson ratio. In view of (9), the above boundary conditions become 

A x + Cj = 0, (11) 

ft^i + ftQ = 0, (12) 

PiCOs(l3ia N )A N + ftsin(fta w )B w + ft cos (/3 2 a N )C N + ft sin (ftfl;v).D/v = 0, (13) 

cos (faa N )A N + sin (P x a N )B N + cos ^a N )C s + sin (fta A ,)D A - = 0. (14) 


As to the continuity conditions between two successive spans, two cases of practical 
interest deserve consideration. 


Case A. Simple support under the rib. In some applications, the flexural rigidity of the 
stiffener is not large enough, and a vertical support is installed under the stiffener to 
suppress the transverse displacement. In this case, the continuity conditions between the 
two typical neighboring spans j and j + 1 are [2] 


or 


w i+ ik +1 =o = 0, 
W;L = „. = 0, 




3 Wj 


dXi 


dw j+1 


dx, 


7+1 


M 


O+Di 


W 

3 3 w.-. 


x j+r 


=0 - M% = (GJ)j i± 1 


dx, +1 dy 


*,+ 1=0 


(15) 


~ D { 


3 2 W: 


7+i 




3 2 w ;+ 1 \ _ / 3 2 W, 3 2 >v,- \ 

+ v ^ + D \ — d. + V i) 

3 y 2 J *,+i=o \ 9 x] 3 y 2 1 


= (GJ), 


3 3 w, 


7+i 


3*,+i3 y 2 


x j+\ 


=0 


where (G/) ; denotes the torsional rigidity of the jth rib. 

Substituting equation (9) into the above conditions of continuity leads to the following 
four equations: 


A j+1 + C /+1 = 0, (16) 

cos (fta ; ).Ay + sin (ft a ; ) + cos (fta ; )C y + sin (fta,)D, = 0, (17) 

-ft sin (ft a,) A, + ftcos(ftfly )ft - ft sin (ftay)C ; + ft cos (fta,)Dy - ftft +1 

-ftDy + i = 0, (18) 

-j8fcos(ftflyMy - ft sin (ft a y ) ft - ft cos (ft a,) C, - ftsin (ftay)ft 

+ ftft +1 + ^ ^ftft + ftc /+1 = ^ftft + i = 0. (19) 
Db 2 D b 2 

Case B. No support under the rib. In this case, the bending, in addition to the torsion, of 
ribs should be taken into account. The conditions of continuity between two consecutive 
spans j and j + 1 read 
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W, 




Bwj 

3x, 


M 


O+i) | 


* 7+1 1 1,*, =0’ 

9 * 7 +i 


dXj+i 

■< 7 + 1=0 

(GJ) r 

9 3 * 7 +i 


or 



3 2 w \ 
+ v i 


/ 3 S+ 1 +i; 3 2 Wy +1 \ 

\3x 2 

IN 

CD 

x r a i 

V 3 xy +1 By 2 I 


= ( GJ) j 5 

*i+ 1=° 3x l+) 3y 




*,+ 1=0 


( 20 ) 


V ( J + \ l+1 =o ~ V% rai = (EI)j 


3 4 w 


i+i 


3v 4 


or 
3 J vv 


dx 


P + (2 — »)■ 


3 2 h>, 


dxfiy 7 


!1sii + ( 2-„)J2vl 

L 3*y +1 


'*,+ 1=0 


= (£/)/ 


3 4 iv 


;+i 


3y 4 


* + 1=0 


3*y +1 3y 2 J 

where lft ;> and F* +1) are the shearing forces in the jth and (/ + l)th spans of the plate; 
( El ) y is the flexural rigidity of the yth rib. 

These conditions of continuity can, in turn, be expressed by the following equations in 
terms of the constants of integration: 

cos (fta ; )Ay + sin (ft a y ) By + cos (ftfly)Cy + sin (ft fly) B, - A J+1 - C j+l = 0, (21) 

~ ft sin (ftfly)Ay + ft COS (ft fly) By - ftsin(ftfly)Cy + ft COS (ft fly) Z), - ftBy + i 

- ftB, +1 = 0, (22) 

-ft COS (ft fly) Ay - ft sin (ftfly)By - ftcos (ftfly)Cy - ft sin (ft fly) By 

+ ftAy +1 + ^ftBy +1 + ftCy +1 + ^ ^ftBy +1 = 0, (23) 

Db 2 D b 2 

ft sin (ftfly)Ay - ftcos(ftfly)By + ftsin (ftfly)Cy - ft cos (ft a,) Z), 

- ^i(|-) 4 Ay +1 + ftBy +1 - ^(^) 4 C J+1 + ftZ)y +1 = 0. (24) 

Introducing the following non-dimensional quantities 


A = 


Pb 2 

v?d' 


r i . ’ 




hZ> 


&B 


N), 


A - Vi 1 + Vi-4 Vi 1 - Vi- 2 ) 


( 25 ) 


the boundary conditions for the simply supported continuous plate equations (11)— (14), can 
be written as 

Aj + C, = 0, (26) 

ftAj + ftQ = 0, (27) 

ft cos (ft r A ,7T)A A , + ft sin (ft r^ir) B N + ft cos (ft wO C 2 + ft sin (ft r^jr) B N = 0, (28) 

cos (ftr N jr)A w + sin (ft r^B^ + cos (ft r^CT, 4 - sin (ft r N jr)D N = 0. (29) 
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The conditions of continuity for Case A, equations (16)— (19). are transformed into the 
following equations: 

A l+l + ft +I = 0, (30) 

cos (ftr ; jr)A ; + sin (ftr ; :r)ft + cos(ftryjr)ft + sin (ft ryff) ft = 0. (31) 

-ft sin (ftryjr)Ay + ft cos(ftr ; jr)ft - ft sin (ft r^) ft 

+ ft cos (ftryff)ft - ftft +1 - ftft* i = 0- (32) 
—ft cos(ftr ; 7r)ft - ft sin (ftrytr)ft - ftcos(ftr y 7r)ft - ft sin (ft r,^) ft 

+ ftA /+1 + r ; jrftJ5 /+ i + ftft+i + T,jrftft +1 = 0. (33) 
For Case B, equations (21)-(24) are rendered into the following form: 

cos (fi l r l TT)A l + sin (ft r,?;) ft + cos (ftryjr)ft + sin (ftr ; jr)ft - A /+1 - ft +1 = 0. (34) 
—ft sin(ftr ; 7r)A ; + ft cos(ftr,7r)ft - ft sin (ft r, ft) ft 

+ ft cos (ftryjr)ft - ftft+i - ftft+i = 0, (35) 
-ft cos (ftrytr)Ay - ft sin (ft r^) ft - ftcos(ftr ; jr)ft - ftsin (ftryjr)ft 

+ ftA y+1 + r ; 3rftft +1 + ftft+i + ryjrftft+i = 0, (36) 
ftsin (ftr ; 7r)A y - ft cos (ft ryff) ft + ft sin (ftr ; jr)ft - ftcos(ftr ; jr)ft 

- £Uy7tAy + , 4- ftft +1 - (UjirCj+i + ftDy+i = 0. (37) 

For a general /V-span continuous plate, we have four equations for boundary conditions 
in the form of equations (26)-(29). Besides, 4 X (A - 1) equations (four equations for 
each rib or interior support such as equations (30)-(33) or equations (34)-(37)) can be 
established from the continuity considerations. Altogether there are 4 x A algebraic 
equations for the same number of unknown coefficients Ay, ft, ft and Dy(; = 1 ~ A) 

[F(A)] 4xW {A} Wx i = 0, (38) 

where elements of matrix [F(A)] are composed of such parameters as those denoted in (25) 
and {A} is a column containing Ay, ft, ft and ft. These equations are linear and 
homogeneous. A non-trivial solution is obtained by setting the determinant of F(A) equal 
to zero, which yields a transcendental equation whose smallest root is the critical buckling 
load A. Having known the buckling load A, equation (38) is used to determine, to an 
arbitrary constant multiple, the coefficients Ay, ft, ft and ft which can then be substituted 
back into equation (9) to obtain the buckling mode shape of the entire plate. Note that, in 
the special case of plates with equally spaced stiffeners, the finite difference calculus 
discussed by Wah and Calcote [5] can be used. In this investigation, however, since we will 
concentrate on the two- or three-span plates with stiffeners not necessarily uniformly 
spaced, the use of the above-mentioned method is not viable here. 


3. A PLATE WITH A SINGLE RIB-STIFFENER 

In order to investigate the variation of the buckling mode of the stiffened plate due to a 
small structural irregularity, here we study the simplest case where there is a single stiffener 
which is slightly misplaced from the mid-span (Fig. 2). This is also the case where the 
discreteness of the stiffener has the most pronounced effect on the buckling load. 
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Let us consider a square plate. Intuitively, we know that the single stiffener should be 
placed as close as possible to the mid cross-section of the plate to produce the highest 
reinforcement on the plate. We will use the following non-dimensional notations for 
specifying the positions of the stiffeners 

r, = i-6, r 2 = i - 6, 6 = *-, (39) 

2 2 a 

where d denotes the misplacement of the stiffener, and 5 is its non-dimensional counter- 
part. When d (or 6) is positive, the stiffener is shifted to the right of its designed position; 
when d (or 6) is negative, the stiffener is located to the left of its designed position. F( A) 
in equation (38) is now an 8 x 8 matrix with elements as follows: 

Case A: F n = 1, F n = 1, F 2 1 = ft, F& = ft, F 35 = ft cos(ftr 2 7r), 

he = ft sin (ft r 2 rr) , F 37 = ft cos (ftr 2 7r), F 38 = ft sin (ft r 2 n) , 

F 45 = cos(ftr 2 jr), F^ = sin(ftr 2 jr), F„ = cos(ftr 2 jr), 

F^ = sin(ftr 2 7r), F 35 = 1, F 57 = 1, F 61 = costftrjjr), 

hi = sin (ft r j7r), F 63 = cos(ftr 1 7r), F M = sin(ftr,7r), (40) 

F 71 = “A sin (A rjjr), F 72 = A cos (ft r, 77 ), F 73 = ftcos(ftr 1 7r), 

F 74 = ft cos (ft /^ tt), F 76 = -ft, F 7g = -ft, F 81 = —ft cos (ft r^), 

F 82 = -ft sin (ftrjjr), F 83 = ftcostftrjjr), F M = -ft sin (ft ^ ft, 

^85 = ft> F^ = Tjffft, Fg/ = ft, Fss = T !jrft. 

The rest of the elements are zero. 

Case B : F n = 1, F 13 = 1, F 2 i — ft, F B = ft, F 35 = ft cos (ftr 2 jr), 

F 36 = ft sin (ftr 2 7r), F 37 = ftcos(ftr 2 jr), F 38 = ftsin(ftr 2 jr), 

F 45 = cos(ftr 2 7r), F^ = sin(ftr 2 jr), F 47 = cos(ftr 2 jr), 

F 48 = sin (ftr 2 J7), F 51 = cos (ftrjff), F 52 = sin(ftr 1 ir), 

F 53 = cos (ft r^), F 54 = sin(ft 2 r 1 7 r), F 55 = -1, F 57 = -1, 

F 6 i = — ftcos(ftrijr), F 62 = ft sin (ft r^), F 63 = -ft cos (ft r^), (41) 

Fm = ftsin(ftrjff), F^ = —ft, F^ = —ft, F 71 = — ftcos(ftr 1 ir), 

F 72 = -ft sin (ft ^jt), F 73 = -ftcos(ftr,jr), F 74 = -ft sin (ftr^), 

F 75 = ~ft, F 76 = ti 7rft, F 77 = ft, F 78 = Tjffft, F 8 ] = ft sin (ft r ]7r), 

F 8 2 — ft cos (ft rj 7r) , F 83 = ft sin (ft r, 77), F*, = -ft cos (ft r, ^), 

F 35 — — (UiTT, F^ = ft, F 87 = —aw, F^ = ft. 

The rest of the elements are zero. The column {A} is now 

{A} = ft, C lf ft, ^ 2 , ft, C 2 , ft} T . (42) 

Setting the determinant of F(A) to zero and using the quasi-Newton root searching 

method, one can find the smallest root of det F(A) = 0, which corresponds to the buckling 
load of the structure. Then, substituting A back into equation (38) and taking any seven 
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Fig. 2. Uni-axial compressed rectangular plate stiffened with a single misplaced rib. 


equations out of the eight equations (38), we can solve for {A}. The buckling mode reads 
wi(*i) = [-4 1 cos (/?!*!) + Bj sin(A^i) 


+ Q cos (P 2 xi) + Di sin (/3 2 *i)] sin 


-(?)• 




w 2 (x 2 ) = [ J 4 2 cos(/3 1 jr 2 ) + B 1 sm(f3 l x 2 ) + C 2 cos (/S 2 jc 2 ) 


+ £> 2 sin(/3 2 x 2 )]sin| 


0 « jc, « — - rf. 
' 2 


4. A THREE-SPAN STIFFENED PLATE 

Now we consider a three-span continuous plate with stiffeners attached to the same 
positions as the interior supports. The plate is all-round simply supported and subjected to 
the uni-axial uniform compression in the direction perpendicular to the stiffeners (Fig. 3). 
For this specific problem, a set of twelve algebraic equations can be established in the form 
of equation (38). The conditions of continuity of the type discussed in Case A are adopted 
here. 

Suppose that the two stiffeners are located at distances fj and § 2 from the left edge, 



Fig. 3. Uni-axial compressed three-span stiffened plate. 
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respectively. Following the same procedure as discussed in the last section, the buckling 
mode for the three spans are expressed as 


w i(*i) = Mi cos (/Si*!) + SjsinC/SiXj) + C 2 c os(/8 3 *i) 
+ Dj sin (&*!)] sin 

^2(^2) = M2 cos (P1X2) + B 2 sin {f$ix 2 ) + C 2 cos (/? 2 x 2 ) 
+ D 2 sin (/J 2 *2)] sin |~"j’ 

w 3 (x 3 ) = [A 3 cos (P\X 3 ) + B 3 sin(y 3 2 x 3 ) + C 3 cos (fi 2 x 3 ) 
+ D 3 sin (/ 3 2 x 3 )] sin 


0 « « 5i,» 

0 =s a ; =s (44) 

0 « *3 « a - £ 2 . 


We are interested in the variation of the buckling load and the buckling mode with the 
small misplacement of the stiffeners. In addition, the optimal position of the stiffeners 
which yields the highest buckling strength is also sought. 


5. NUMERICAL RESULTS AND DISCUSSION 

Numerical calculations are performed for both the single rib-stiffened plate and the 
stiffened three-span continuous plate. Structures with different parameters for the torsional 
and flexural rigidities are also investigated. 

For the plate with a single stiffener, attachment of the stiffener to the mid cross-section 
of the plate provides the structure with the most favorable load carrying capacity. This 
conclusion holds true for both Case A and Case B. It is found that for Case A, where 
there is a support which prevents the vertical displacement of the plate, the magnitude of 
the nondimensional torsional rigidity r has only a moderate effect on the buckling load 
when r is larger than 10 (Fig. 4). Deviation of the stiffener from its supposed mid-span 
position reduces the buckling strength, but more importantly, it changes the buckling mode 
from an overall buckling of the plate to a local buckling of the plate segment with longer 
span. The more misplaced the stiffener, the greater the reduction in buckling load will be, 
and the more localized the buckling mode becomes, that is, the deflection of the plate on 
one side of the stiffener is much greater than that on the other side. For example, for Case 
A with a torsional rigidity of r = 20.0, the ratio of the maximum deflection in the left 
segment to that in the right segment is 4.5 when <5 = 0.01. If a bigger misplacement is 
involved, say <5 = 0.02, then the ratio of the maximum deflections in the two segments 
increases to 7.0. However, for the stiffener with torsional rigidity t < 5, small misplacement 
does not significantly affect the buckling load; for instance, when r = 2, a deviation of 
magnitude <5 = 0.05 produces only 4% reduction in buckling load. Figures 5 and 6 show the 
buckling mode shape of the plate in Case A for different values of r. With a stiffener of r 
larger than 30, the shorter segment of the plate is almost undeflected as buckling mode is 
localized in the longer segment. For Case B, the flexural rigidity of the stiffener plays a 
more important role in the buckling strength than the torsional rigidity, although the 
influence of torsional rigidity is still remarkable in the buckling mode shape. For example, 
it can be seen from Fig. 7 that only the stiffener with flexural rigidity a> ^ 5 has noticeable 
strengthening effect, and that when c 0 falls below two, the position of the stiffener becomes 
almost irrelevant for the magnitude of the buckling load. When a plate is reinforced with a 
rib of moderate flexural stiffness, the longer segment of the plate is severely deflected at 
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Rib misplacement parameter 6 


Fig. 4. Loci of buckling loads for a plate stiffened by a single rib with different values of r (Case A). 



Fig. 5. Buckling mode localization for a plate stiffened with a single rib of r = 20.0 (Case A). 



Fig. 6. Buckling mode shape for a plate stiffened by a single rib with misplacement 6 = 0.02 (Case A). 
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-0.10 -0.06 -0.02 0.0Z 0.06 0.10 

Rib misplacement parameter 6 

Fig. 7. Loci of buckling loads for a plate stiffened by a single rib with different values of to (r = 10.0, Case B). 

the onset of buckling while the short segment experiences only a scant deformation. So the 
buckling is still fairly localized, as can be seen from the mode shapes depicted in Figs 8 and 
9. It is interesting to note that the cut-down in the overall strength of the plate by 
mispositioning a stronger stiffener can be greater. For instance, when r = 30 and to = 20, a 
5% deviation from the mid-point produces 13% decrease in the buckling strength. Thus we 
can see, a unilateral increase in the stiffener’s strength may make the whole structure 
highly sensitive to the misplacement (which can be regarded as a special kind of initial 
imperfection) in the sense that a small misplacement of the stiffener or interior support can 
lower the buckling load of the plate, and more importantly, localize the buckling mode 
shape. 

For the three-span continuous plate, numerical results show that, as far as the buckling 
load is concerned, equally spaced stiffeners such that the three spans of the plate have the 
same length are not most beneficial. As compared with the single rib-stiffened plate, the 
three-span plate is even more sensitive to the misplacement of the stiffener. For example, if 
one stiffener is fixed at ^ = a/3, the other stiffener is supposed to be located at § 2 = 2a/3 
but somehow it deviates from this position by, say, ^ = -0.02 (the negative sign 
representing the misplacement is in the negative x direction), the buckling load is 



Fig. 8. Buckling mode localization for a plate stiffened by a single rib of x = 20.0 and to = 10.0 (Case B). 
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decreased by 9.5% from its counterpart without misplacement. Interestingly enough, some 
patterns of the misplacement are detrimental, while others can be helpful. For instance, 
suppose the stiffeners are designed to be located at ci = a/3 and g 2 = 2a/3, respectively. 
The combination of the misplacement by the magnitude = & = 0.02 (meaning the 
misplacement are all in the positive x direction) from = a/3 and g 2 = 2a/3, respectively, 
cuts down the buckling load by 9.6%. However, misplacement of dj = -& = — 0.02 
(meaning the left stiffener has been moved slightly to the left and the right stiffener to the 
right) boosts the buckling strength by 11% (Fig. 10). As for the buckling mode, in the 
majority of the situations, the buckling is still severely localized (Fig. 11). This is shown by 
Fig. 12 where a small misplacement = —0.01 of the right stiffener triggers the onset of 
the local buckling in the third span of the plate at a lower load than its counterpart without 
the misplacement. When this happens, the other two spans hardly deflect at all. 

Numerical analysis shows the optimal stiffener layout for stiffeners with torsional rigidity 
r = 20 is £i = 0.329a and § 2 - 0.671a. The corresponding buckling load is 19.6% above the 
buckling load with two identical stiffeners positioned at = a/3 and g 2 = 2a/3. This result 
can also be interpreted as follows: for misplacement dj = 0.334 — 0.329 = 0.005 and 
62 = 0.667 — 0.671 = -0.004, the buckling load is decreased by 16% ((1 — 1/1.19) x 



Fig. 9. Buckling mode shape for a plate stiffened by a single rib with misplacement <5 = 0.01 (Case B). 



Fig. 10. Effect of misplacements on the buckling load of a stiffened, three-span plate (ri = t 2 = 20.0). 
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Fig. 11. Buckling mode shape for a stiffened, three-span plate without stiffener misplacement (r, = r ; = 20.0. 

b = 0 . 0 ). 




Fig. 12. Buckling mode shape for a stiffened, three-span plate with slight stiffener misplacement (ft — r 2 — 20.0, 

6 , = 0 . 0 , &2 = - 0 . 01 ). 

100%). This again demonstrates the high sensitivity of optimally designed structures to 
small imperfections. This phenomenon was discussed by Budiansky and Hutchinson [18] as 
well as Zyczkowski and Gajewski [19]. Moreover, it is found here that the^ optimal pattern 
of the stiffener layout is almost independent of the specific value of r, as long as the two 
stiffeners are identical. For example, even when the torsional rigidity r is decreased to 5.0, 
the most favorable positions of two stiffeners are hardly changed, as they are now 
§! = 0.330a and § 2 = 0.670a. Figure 13 depicts the buckling mode for such a situation, 
from which we can observe that with the optimal stiffener layout the buckling mode is a 
global one, that is, all parts of the plate deflect to comparable degree and the potential 
capability of the structure is fully tapped. 

In conclusion, the buckling mode localization phenomenon due to small misplacements 
in the stiffened or continuous plates should not be overlooked, especially in those 
applications where the mode shape is a significant concern. Due to the imprecision in the 
fabrication, the misplacement is always present and can to the large extent affect the 
buckling characteristics of the structure. When the structure is designed in terms of the 


1530 


I. ELISHAKOFF et al. 


CD 

o’ 



Fig. 13. Buckling mode shape for a stiffened three-span plate with optimal stiffener placement (r, = t, = 20.0. 

I, = 0.329a, h = 0.671a). 


optimal stiffener layout, misplacement reduces the buckling load and causes the buckling 
mode to become highly localized in a manner that one segment of the structure deflects 
appreciably while the capability of the other parts of the structure has not been brought 
into full play. 
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Abstract — In this paper we present a closed-form solution for vibrational imperfection sensitiv- 
ity — the effect of small imperfections on the vibrational frequencies of perturbed motion around the 
static equilibrium state of Augusti’s model structure (a rigid link, pinned at one end to a rigid 
foundation and supported at the other by a linear extensional spring that retains its horizontality as 
the system deflects). We also treat a modified version of that model with attendant slightly different 
dynamics. It is demonstrated that the vibrational frequency decreases as the initial imperfections 
increase. 

Keywords: dynamics of rigid links, influence of imperfections on link vibrations, rigid body 
dynamics 
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INTRODUCTION 

The theory of dynamic imperfection sensitivity was put forward in a number of studies by 
Budiansky and Hutchinson [1-3] as a counterpart to the static framework developed by 
Koiter [4], Whereas the latter investigation yielded the maximum (limit) load which can be 
sustained by the structure under static conditions, the former studies resorted to a step-load 
applied to the system in order to find the effect of unavoidable small imperfections on the 
maximum dynamic load. A general framework for the problems of suddenly loaded 
structures was considered by Simitses [5], 

In recent years, dynamic imperfection sensitivity has also been studied in a different 
context, namely, with the system loaded statically up to some specific load level at which it 
is dynamically perturbed. The resulting motion possesses frequencies which are influenced 
both by the load level and by the initial imperfections. (For details, the reader is referred to 
the papers by Rosen and Singer [6], Singer and Prucz [7], Elishakoff et al. [8], Hui and 
Leissa [9, 10], Shilkrut and Virlan [11], Nash [12] and Valishvili [13].) Here, the general 
behavior can be elucidated by examining the various effects involved on model structures. 
One such model structure was pioneered by Budiansky and Hutchinson [1-3], and the 
influence of the initial imperfections and its vibrational frequencies was studied by 
Elishakoff et al. [14], Another interesting model in the static context was thoroughly 
investigated by Augusti [15] and by Thompson and Hunt [16]; Yizhak et al. [17] studied 
the dynamic buckling loads of this model, with relevant closed-form expressions. Other 
models were considered by Souza [18] and by Kounadis et al. [19], We will now deal with 
the original Augusti model and its modified version, with regard to their vibrational 
imperfection sensitivity in the vicinity of a non-linear static state. 


ANALYSIS 

Consider the cantilever shown in Fig. 1(a); it consists of a rigid link (of length L) pinned to 
a rigid foundation and supported by a linear extensional spring (of stiffness k) which retains 
its horizontality as the system deflects. The initial imperfection is modeled by the angular 
deflection a 0 from the vertical position, and the non-dimensional total displacement is 
denoted by x. We first consider the static case. 
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Fig. 1. 


Due to the equilibrium, we have 

£M a = PxL - k(x - x 0 )y/l -x 2 L = kL [ax - (x - x 0 )y/l -x 2 ] = 0 (1) 

where 

<r = , Pci = kL, (2) 

*cl 

P ci being the classical buckling load. From equation (1): 

(3) 

The limit load is given as the maximum load the system can support, so that 

da da x 0 - x 3 

d* ’ dx xV 1 -* 2 

leading to 

x 5 = Xo /3 (5) 

a s = (1 - x 3 ' 3 ) 3 ' 2 (6) 

which coincides with the formulae of Augusti [15] and Thompson and Hunt [16]; x s is the 
total displacement at which the static limit load a s is achieved. 

Let us now consider the dynamic case. At some load level a, which corresponds to the 
total displacement x, we superimpose the perturbation w(t), so that the total displacement 
becomes 


x(t) = X + w(t). 


( 7 ) 



The equation of motion reads 
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- M a = moiL 2 (8) 

where the dot denotes differentiation with respect to time. 

Denoting 

a>o = k/m, (9) 

co 0 being the natural frequency of the mass m under small motion. Equation (8) is rewritten as 

coo 2 a = ox — (x ~ x 0 )-»/l — * 2 (10) 

where (see Fig. 1): 


x = sin a 


a = • 


: + - 


(x) 2 


V'l - X 2 (1 - 2C 2 ) 3/2 
We seek a solution in the form 


= (l-x 2 ) 


- 1/2 


w + 


(w ) 2 

1 -x 2 


( 11 ) 

( 12 ) 


w(t) = W e imt 

subject to restrictions such that instead of equation (12), we arrive at: 


coq 2 w = [ ax — (x — x 0 ),/l — if 2 ] ^/l — x 2 . 
With equation (13) in mind, we have 

O)o 2 w= —Q 2 W e icot 

where 


Q = 


co 

CD 0 ' 


Denoting 


f{x) = (x- Xoj^/l - x 2 , 
equations (1) and (14) are rewritten as follows: 

ax =/(x) 


^ 2W = lf(X ~ ffX] y/ 1 — X 2 . 

We note that at the limit load, and after linearization, we obtain: 

do- /'(x) - a 


0 = 


dx 


or, by equation (3): 


= [/'(*) -^Vi -x 2 




(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

( 20 ) 
(21) 

( 22 ) 


Noting that at the limit load dc/dx = 0, and by inspecting equation (22), we arrive at the 
important conclusion that at that level the natural frequency Q vanishes identically: 

a = a s : Q 2 = 0. (23) 

Bearing in mind equations (4)— (6), we obtain the following interesting representation: 

Q 2 =^ (24) 

which must be solved in conjunction with equation (3). 

The following question now arises: “How does the non-dimensional squared frequency 
Q 2 change, for fixed loading a = const (cr < 1) as the imperfection increases from x 0 = 0 (i.e. 
Q 2 = 1 — a) to x 0 = (1 — a 213 ) 312 , at which O 2 = 0?” 
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From equation (3) we deduce 


Q 2 = 1 - x 2 


CJ 

y/l-x 2 


which is a decreasing function of x. Also, observing that 


and 


we conclude that 


dtr 

dx 

do 

dx 0 


x 0 — x 


x 0 = const X 2 ^/ 1 — X 2 




>0 


<o 


dx 

dx 0 


>0 


(25) 


(26) 

(27) 

(28) 


and 


d(Q 2 ) 

dx 0 


<0 


(29) 


implying that Q 2 decreases monotonically. 

In order to interrelate Q 2 , x 0 and o, we recast equation (3) in the form of the polynomial 
equation: 

p 4 (x) = x 4 — 2x 0 x 3 — yx 2 + 2x 0 x — x§ = 0 (30) 


where 


7 = 1 ~a 2 ~xl. 

(31) 

Now, for a < 1, x 0 < 1, we arrive at 


a 2 + x o < <r 2/3 + x 0 2/3 < 1 

(32) 


that is, y is always positive. 

Now equation (24) also can be written as a polynomial equation: 

q 3 {x) = x 3 + Q 2 x — x 0 = 0. 


(33) 


We are looking for the conditions under which equations (30) and (33) have a common 
root. Equation (30) could be put as or: 

p 3 {x) = — p 4 (x) + xq 3 (x) = 2xqX + (y + Q 2 )x 2 — 3x 0 x + Xo = 0 (34) 

q 2 (x) = p 3 (x) + 2 x 0 q 3 (x) = (y + Q 2 )x 2 — (3 + 2fi 2 )x 0 x + 3 xq = 0. (35) 

Instead of equation (33), we resort to 

q 2 (x) = [? 3 W^o + p 3 (x)]/x = 3x 0 x 2 + (y + Q 2 )x + (Q 2 — 3)x 0 = 0. (36) 

For q 2 (x) = 0 and q 2 (x) = 0 to have a common root, it is necessary and sufficient that 
their resultant R{q 2 ,q 2 ) vanish: 


R(q 2 ,h) = 


y + Q 2 — (3 + 2f2 2 )x 0 
0 y 4- Q 2 

3x 0 y + Q 2 

0 3x 0 


3xq 

- (3 + 2Q 2 ) 
(Q 2 - 3)x 0 
y + 12 2 


0 

3x 0 2 

0 

(Q 2 - 3)x 0 


= 0 


(37) 


which yields the final equation 

2(xq - a 2 ) - 27x^A 2 + y 3 + (y 2 - 3y)Q 4 + Q 8 


0. 


(38) 
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Two particular cases associated with this equation are: 

(a) ' x 0 = 0: Q 2 = 1 — cr (39) 

(b) cr = 0: Q 2 = 1 — Xq . (40) 

Figure 2 shows the dependence of the non-dimensional frequency upon the initial 
imperfection parameter x 0 , for different values of the loading a. For a zero initial imperfec- 
tion, the curves initiate at (1 - cr) and descend as the initial imperfection increases. Figure 3 
shows the non-dimensional frequencies vs the non-dimensional loading, for different initial 
imperfections. 



Fig. 2. 



Fig. 3. 
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MODIFIED AUGUSTI MODEL 

Consider now' a slightly modified Augusti model, statically equivalent to the original 
version but with different dynamics [Fig. 1(b)]. The equation of motion now reads 


mxL = P 




(41) 


where F is the force in the spring. 

This equation could be rewritten as 

a>o 2 x = 


ax 




(x - x 0 ) 


(42) 


where x 0 is again the imperfection. 

Equation (42) in conjunction with (13) yields 

Q 2 = - 


+ 1. 


(43) 


(Vl -x 2 ) 2 

This expression again reveals that Q 2 is a decreasing function of x; equivalent to the latter 
equation is 

,3 


Q 2 = 


= 1 


or 


Q 2 = 1 - 


We make the following substitutions: 


x — x 0 

x(l — x 2 ) x x(l — x 2 ) 
(x — x 0 )x a 


1-x 2 


yr 


(44) 

(45) 



a = Xo /3 = x 0 

(46) 


a a a 

(l-a) 3 ' 2 ’ 

(47) 

Hence 

Q 2 = 1 2 ~ a 

z(l - az 2 ) 

(48) 

yielding 

(1 — Q 2 )az 3 + Q 2 z — a = 0. 

(49) 


On the other hand, from equation (3) we arrive at 

(z — a) 2 (az 2 — 1) + A(1 — a) 3 z 2 — 0 


(50) 


or 

az 4 — 2a 2 z 3 — /? 2 z 2 + 2 az — a 2 = 0 (51) 

where 

ft = l — a 2 — A 2 (l - a) 2 . (52) 

Figure 4 shows the non-dimensional natural frequency squared vs the non-dimensional 
loading a. As a approaches the critical values a s , the non-dimensional frequency vanishes 
identically. At a — 0 all curves initiate at 1 — Xq. Figure 5 shows the behavior of Q 2 as 
a function of x 0 , and here again, as in Fig. 2, the natural frequencies decrease as x 0 increases. 


CONCLUSIONS 

We investigated the effect of initial imperfections in Augusti’s model structure, in the 
vicinity of the non-linear static state. Qualitatively, the effect is analogous to that associated 
with buckling problems, where the buckling load decreases as the initial imperfections 
increase. Accordingly, the concept of vibrational imperfection sensitivity is seen to be very 
important for engineering applications, since the vibrational frequency actually experienced 
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by the system could be much lower than the design level. In this paper, the closed form 
solutions have been derived to serve as benchmark solutions for the vibrational imperfec- 
tion sensitivity. 


Acknowledgement — Support provided by the NASA Langley Research Center through grant No. NAG-1-1310 is 
gratefully acknowledged and appreciated. 


REFERENCES 

1. B. Budiansky and J. W. Hutchinson, Dynamic buckling of imperfection sensitive structures. In Proceedings 
of the Eleventh International Congress of Applied Mechanics, H. Gobrtler (ed.), Munich, Germany, 
pp. 636-651 (1964). 




236 


I. Elishakoff et al. 


2. B. Budiansky, Dynamic buckling of elastic structures: criteria and estimates. In Dynamic Stability of Structures, 
G. Herrmann (ed.), pp. 83—106. Pergamon Press, Oxford (1967). 

3. J. W. Hutchinson and B. Budiansky, Dynamic buckling estimates. A1AA J. 4 , 525-530 (1966). 

4. W. T. Koiter, On the stability of elastic equilibrium, Ph.D. thesis. Delft University of Technology, N.J. Paris, 
Amsterdam (in Dutch). English translations: (a) NASA TTF-10, 833 (1967); (b) AFFDL-TR-70-25 (1970). 

5. G. J. Simitses, Dynamic Stability of Suddenly Loaded Structures. Springer, Berlin (1990). 

6. A. Rosen and J. Singer, Effect of axisymmetric imperfections on the vibrations of cylindrical shells under axial 
compression. AlAA J. 12 , 995-997 (1974). 

7. J. Singer and J. Prucz, Influence of initial geometrical imperfections on vibrations of axially compressed 
stiffened cylindrical shells. J. Sound Vibr. 80 , 117-143 (1982). 

8. I. Elishakoff, V. Birman and J. Singer, Small vibration of an imperfect panel in the vicinity of a nonlinear static 
state. J. Sound Vibr. 114 , 57-63 (1987). 

9. D. Hui and A. Leissa, Effects of uni-directional geometric imperfections on vibrations of pressurized shallow 
shells, lnt. J. N on-Linear Mech. 18 , 279-285 (1983). 

10. D. Hui and A. Leissa, Effects of geometric imperfections on vibrations of biaxially compressed rectangular flat 
plates. J. Appl. Mech. 50 , 750-756 (1984). 

11. D. I. Shilkrut and P. M. Virlan, Dynamic study of the stability of all equilibria of geometrically nonlinear 
shallow spherical shells. Stroitel'naya Mekhanika i Raschet Sooruzhenii, pp. 54-58 (1975) (in Russian). 

12. W. A. Nash, Free vibrations of initially imperfect cylindrical shells. In Proceedings of the 8th Canadian 
Congress of Applied Mechanics, Moncton, pp. 365-366 (1988). 

13. N. V. Valishvili, Computer Aided Design of Shells of Revolution, pp. 265-271. Mashinostroenie, Moscow (1976) 
(in Russian). 

14. I. Elishakoff, V. Birman and J. Singer, Effect of imperfections on the vibrations of loaded structures. J. Appl. 
Mech. 51, 191-193 (1984). 

15. G. Augusti, Stability di Strutture Elastiche Elementari in Presenza di Grandi Spostamenti, Publication No. 
172, Department of Structural Sciences, Faculty of Engineering, University of Naples, Italy (1964). 

16. J. M. T. Thompson and G. W. Hunt, A General Theory of Elastic Stability. John Wiley, London (1973). 

17. E. Yithak, I. Elishakoff and M. Baruch, Dynamic imperfection sensitivity of a simple spring model structure. 
J. Mech. Struct. Machines 16 , 187-199 (1988). 

18. M. A. de Souza, Dynamic behaviour of structural elements liable to buckling. In Proceedings of the IUTAM 
Symposium on the Inelastic Behaviour of Plates and Shells, PUC/RJ-Rio de Janeiro, Brazil, 5-9 August (1985). 

19. A. N. Kounadis, .J. Raftoyiannis and J. Mallis, Dynamic buckling of an arch model under impact loading. 
J. Sound Vibr. 134 , 193-202 (1989). 



Pergamon 



Chaos, Solitons & Fractals Vol. 7, No. 8, pp. 1179-1186, 1996 
Copyright © 1996 Elsevier Science Ltd 
Printed in Great Britain. All rights reserved 
0960-0779/96 $15.00 + 0.00 


0960-0779(95)0110-7 


Imperfection Sensitivity Due to the Elastic Moduli in the 


Roorda-Koiter Frame 




I. ELISHAKOFF and Y. W. LI 

Department of Mechanical Engineering, Florida Atlantic University, FL 33431-0991, USA 


and 


J. H. STARNES JR 


->S -.3 7 


J/o 8 9 'S 


NASA Langley Research Center, Hampton, VA 23665-5225, USA 
(Accepted 13 November 1996) 


/ f 


Abstract— In this study, it is demonstrated through a simple example of the Roorda-Koiter frame 
that the unavoidable dissimilarity in the distribution of elastic moduli may further reduce the 
load-carrying capacity in addition to the well-recognized effect of initial geometric imperfections. 
Copyright © 1996 Elsevier Science Ltd. 


INTRODUCTION 

The works by Koiter [1, 2], Budiansky and Hutch i nson [3] established that the initial 
geometric imperfection— deviation from the nominally ideal shape— plays a major role in 
the reduction of the load-carrying capacity of structures. An additional source of such a 
reduction was identified recently by Koiter et al. [4, 5] as a non-uniform thickness of the 
structure. To the best of the authors’ knowledge, there is only a single work, by Ikeda and 
Murota [6], that points out that the symmetry breaking in the distribution of elastic moduli 
may cause additional reduction in the loads the system is able to sustain. In particular, they 
analyzed the von Mises truss in such a context. In this note, we perform analytical and 
numerical investigation on a frame structure to study the effect of dissimilarity of elastic 
moduli on the load-carrying capacity of the structure. 


ANALYSIS 

In 1965, Roorda [7] conducted a set of experiments on the two-bar frame, subjected to 
an eccentric load. Koiter [8] performed an initial-imperfection analysis of this frame, 
referred to in the literature since then as the Roorda-Koiter frame [9]. Roorda and Chilver 
[10] and Bazant and Cedolin [11] analyzed this structure from different perspectives; El 
Naschie studied a discrete analog of the Roorda-Koiter frame [12] (for additional 
references, see also the text of Brush and Almroth [13]). In this note, the vertical and the 
horizontal segments of the two-bar frame (Fig. 1) may have different lengths L x and L 2 
and different Young’s moduli E x and E 2 - Our aim is to study the effect of a small 
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difference in the elastic moduli between the two segments, namely, the vertical and the 
horizontal bars, of the frame structure. It is such a difference that may trigger the 
dissimilarity in the elastic moduli of the structure. 

Following those previous investigators, we first analyze the perfect structure. 

The total potential energy of the two-bar frame is 


n 


f 

Jq 


+ 


i e x a 

[ d § 

l 2 

. dr 


iwie-en - 


JO f 2 

The change in potential energy is obtained by letting 

£ lo + f> V Vo + »?, v -* v 0 + v. 


w — > w 0 + w 


( 1 ) 


( 2 ) 


where 


% = v 0 = w 0 = 0, 



(3) 


represents an equilibrium state on the primary equilibrium path in the neighborhood of the 
bifurcation point A. The energy expression may take the following form 


where 



(4) 


( 5 ) 



Elastic moduli 
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1 dv y 
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r/ 2 

J 2 dy ' 

vd7/ 


1 o4 


4! 


<5 4 n 


r^Ej^A ( dwY , E 2 A f L *l dti \ 4 , 

— - 1 dx + — dy 

Jo 8 \ dx / 8 Jo \ dy J 


( 6 ) 

(7) 


correspond to second, third and fourth variations of the potential energy, respectively. 

According to Koiter’s initial postbuckling theory [1], in the initial postbuckling range, the 
incremental displacement components are of the form of the classical buckling mode 

I = Pa% 1, n = PaVu V = Pa»u w = PaW 1 (8) 

where £ l5 rj 1 , v l , w t are normahzed classical buckling load. fi A is the rotation angle at 
bifurcation point A, and could be viewed as the amplitude parameter. 

Substituting the above into the second variation leads to 


i# n = y ^(— y + ^y 

2! Jo \ dx / \ dx ) 2 \ dx 2 ) 

Jo \ dy / \ dy 2 / J 


dxp\ 


(9) 


Performing variation on the above expression, we obtain, according to the Trefftz 
criterion, 


= [ l \e a (— - p( dWl \ d ( gw, i) + ElI ( d2w i \ d2 (^i 
\2! / Jo . 1 \dx/ dx \ dx J dx 2 \ dx 2 / dx 2 


d XpA 


+ 


+b4—J 

Jo L \ dy / dy \ dy 2 / 


d 2 (<5ni) 


dy / dy 

Integration by parts and rearrangement gives 

|3 


dy 2 / dy 


(10) 


d ypA = 0. 


La— 5§x 

L 

+ 

L 2 A^6ijx 

'' - [( 

dx 

0 

dy J 

0 L\ 


. d w, 

EJ - + P 


E 2 I^l8v 1 

dy 3 


+ 


E l I 


dx 


d 2 Wi d(6vy 1 ) 


d "'> 
dx 


dx 2 dx 

l 

dy 2 


■ Lj 

+ E 2 I 
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d 2 vi d(dn 1 ) 


dy 2 dy 


_ {\ A iE dil i X - 

Jo dx 2 Jo dv 2 Jo l 4c 4 dx 2 I 

L 


(ii) 


l 2 d 4 n, 

+ | E 2 I -dv 1 dy = 0. 


dy 


The satisfaction of the above equation results in boundary conditions (using £|* =Ll = 

-H 


b = L 2 ’ V\y=L2 w \x=L 2 )’ 

d 2 w l 


dx 2 


= 0 at x = 0; 


d 2 U! 

dy 2 


= 0 at y = 0 


,d|i 


E,A-^ + E 2 I 
dx dy 


d 3 iq d 2 w 1 d 2 Uj 


dx 2 


dy 2 


= 0 at x = L 1; y 


E ? A 


drji 
dy 

and governing equations: 


EJ 


d 3 Ux 

dx 3 


,dw\ 


- P = 0; at x = Lx, y 

dx 


( 12 ) 

(13) 

(14) 
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d 2 £i 
dx 2 
d 2 ??i 
dy 2 

Integration of the above equations gives 

A»fi = Q- " 


= 0, 

E X I- 

= 0, 

d 4 U! 

^,,4 


d 4 wq p _cPw x _ 


dx 4 
= 0. 


dx 2 


PaVi - C 2 


E X A E 2 A 

Pa^i = C 3 sin(kx) + C 4 x, p A v x = C 5 y + C 6 y 3 

where k 2 = P/E X I, and C, (i = 1, 2, . . 6) are integration constants. 

Using boundary conditions at x = L 1; y — L 2 : 

^2 „ 


w = rj: C 3 sin(A:L 1 ) + C 4 L X 


§ — n: C5L2 + C 6 L 2 — 


£2^ 


£ x .4 


C x 


dw 

dx 


= — p A \ kC 3 cos(kL x ) + C 4 


'Pa 


dv 

d y 

,d?i 


— —p A : C 5 + 3 c 6 L 2 — —ft/ 
r d 3 ui 


£ X A— + £ 2 7 — -i = 0: C x + 6 E 2 IC 6 = 0 
dx dy 3 

— + — — = 0: — & 2 C 3 sin(£L x ) + 6C 6 L 2 = 0 

dx 2 dy 2 


E 0 A 


dr? x 


Enl 


d 3 u x 


,dvvi_ 

dx 


0: C 2 - P[kC 3 cos (kL x ) + C 4 ] = 0. 


(15) 


(16) 

(17) 

(18) 

(19) 

( 20 ) 
( 21 ) 
( 22 ) 
(23) 


dy “ dx 3 

The above seven equations are linear, homogeneous algebraic equations in terms of j} A , 
Ci, . . ., C 6 . The non- triviality condition leads to 


6(£L X ) cos (£L X ) 


(£L x ) 2 sin(£L x ) — 6sin(£L x ) - 2(£L x ) 2 -^-sin(£L x ) 

e x al\ l x 


+ 


E X A E X AL X L 2 


- -- (kL x ) 2 sin (kL x ) + 2(kL x ) 2 P 


E 2 a 




sin(/c£ x ) = 0. (24) 


Since our discussion is limited to the elastic range, terms containing P/E X A, 
6E 2 l/E x ALl, P/E 2 A and 6 E x I/E x AL x L 2 are negligibly small, and could be omitted in the 
first approximation. Thus, we have an approximate characteristic equation as the following 


tan(k£ x ) = ^ (25) 

1 + \{kL x ?±E 

E 1 

from which, the critical buckling load parameter k d can be determined. The normalized 
classical buckling mode is found to be 
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m 


W x = 


3 Pc 

E\A(k d Li ) 2 

3 P c i 

E l A{k cl L 1 f ] 
3 


fi 


(M-l ) 2 

= b( 1 - 


x — Li 

y 2 \ 


sin (/c d x) 


(26) 


sin(A: c ;L 1 
P = P d + (A - 1)P C/ . 


For the case L l = L 2 = L, the buckling load parameter is found to be k d L = 3.72. 
Substituting the above expressions, we have 


— &TI = (o. 478 P c/ L + 0.109— P C ,L - 0.587A]^. (27) 

2! \ E x ) 

If we assume that there is a dissimilarity in the distribution of elastic moduli, i.e. 
generally E x + E 2 



E 2 = E 1 ( 1 + s) 

(28) 

then equation (29) reads 

— = 0.587(1 + 0.186s - A )P cl Lp 2 A . 

(29) 


2! 


The third variation is 

^-<5 3 II = 0.149 P d Lp A . 

(30) 


The approximate expression for the potential energy increment All is the sum of the 
second and third variations of the energy expressions. For equilibrium, d(An)/d/3 /1 = 0. 
For p A =£ 0, we obtain 


A = 1 + 0.186s + 0.381 /? a 


(31) 


or, we can re-write equation (31) as 

K = — = 1 + 0-381 Pa-, P'd = PM + 0.186s). (32) 

P’ cl 1 + 0.186s 

We now proceed with the analysis of the geometrically imperfect structure. For an 
eccentric load applied at a distance (j)L to the right of point A, the potential energy 
expression must be modified by adding a term [13] 

= ~PL<p/3 A = -X x P' cl L<t>f} A (33) 

and the second variation and the third variation take the form 


= a 2 (i - a dP'Mfa 

^-<5 3 n = A 3 (l - A )P' cl Lp 3 A 

where A 2 = 0.587/(1 + 0.186e) 2 and A 3 = 0.149/(1 + 0.186s). Then 
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An = itfn + i-6 3 n + q 0 = [(1 - A } )A 2 p 2 A + A$\ - iMAP'clL. 

For the equilibrium, we have &(&Xl)/&p A = 0: 

3A 3 


Aj = 1 + A — -^-A, 

2A 2 2A 2 p A 


or, in another form, 


KPa = Pa + ^P\ 


-<PK 


(35) 


(36) 


(37) 


2A 2 2^4.2 

Differentiating both sides of the above equation with respect to , we have 
a dAj , 3A 3 1 dA x 

+ Aj - 1 + -—Pa ~ — • (38) 

d/?A ^2 2y4 2 d/3^ 

Setting dX l /&p A equal to zero, we obtain the limit-point load factor k L for the imperfect 
structure 


A L = 1 + ^p A 


o _ (2l ~ 1)^2 

* 3A 3 

Substituting back into equation (36), we have 

A l = 1 - (3A 3^) 1/2 ( - 0)1 / 2 , 

^2 

As a first approximation, we can set X L = 1 at the right-hand side of the above equation, 
and we arrive at 

A L - 1 - (3A3 -- 2 - (-0) 1/2 = 1 - «(-tf>) 1/2 , a = (42) 

-4-2 -A 2 

where a is a parameter describing the imperfection sensitivity of the structure. For e = 0, 
a = 1.15; e = 0.1, a = 1.17; £ = —0.1, a = 1.11. Using these data, we can plot the 
relationship between the buckling load and the initial imperfection parameter. This is 
depicted in Fig. 2, where the three curves correspond to three different cases of elastic 
moduli. 

For the case where L 2 /Lj = 0.5, the buckling load parameter obtained from equation 
(25) is ( k d L x ) = 3.97. Following the same procedure as above, we again end up with 
expression (42). However, A 2 and A 3 now have the following expressions 


(39) 

(40) 

(41) 


^2 


0.37 

(1 + 0.5l£) 2 ’ 


At, 


0.09 

1 + 0.5l£ 


(43) 


For £ = 0, a = 1.40; e = 0.1, a = 1.51; s = —0.1, a = 1.30. Having these, we can plot the 
curves of buckling load reduction as in Fig. 3. 


DISCUSSION 

Figures 2 and 3 demonstrate the usual pattern of the buckling load reduction due to 
initial imperfection. However, what is worth noticing here is a further reduction in buckling 
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Fig. 2. Buckling load reduction {L 2 /Li = 1.0). 



Fig. 3. Buckling load reduction (L 2 /Li = 0.5). 


load because of the presence of dissimilarity in the elastic moduli and the geometric 
configuration. When e is positive, it means that the horizontal bar possesses a bigger elastic 
modulus than the vertical bar, which results in a stiffer structure. Likewise, when L 2 /Lj is 
less than unity (the horizontal bar is shorter), the overall stiffness of the structure is again 
increased. On the one hand, a stiffer structure has a larger buckling load. On the other 
hand, Figs 2 and 3 show that the stiffer the structure is, the more sensitive it is to the initial' 
imperfection. As we can see also from those two figures, when the stiffness is decreased 
either through a reduced elastic modulus (e< 0) or from an increased length (L 2 /L x > 1), 
the structure becomes less sensitive to the initial imperfection. The change in the overall 
stiffness of the structure comes, as is shown here, from the dissimilarity. Thus, one may 
conclude that the dissimilarity in elastic moduli could contribute to the change in sensitivity 
of the structure to the initial imperfection. One may intentionally introduce such a 
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dissimilarity in elastic moduli so that the structure has a decreased sensitivity to initial 
imperfections. In this simple two-bar structure, the decrease in sensitivity due to dissimilar- 
ity in elastic modulus is not remarkable when the two segments of the frame have the same 
length, namely, L x = L 2 . However, for other geometric configurations, the effect may 
become more pronounced. For L 1 = 2 L 2 (Fig. 3), say, at <j)= 0.25, A = 0.3 for e = 0, and 
A = 0.24 for e = 0.1, which indicates 20% decrease in the limit load; for e = —0.1, A = 0.3, 
which amounts to 17% increase in load-carrying capacity. It appears that the subject is 
worth pursuing in the direction of shell structures to see how important the effect of 
non-homogeneity of elastic moduli is on the imperfection sensitivity. 

Acknowledgement — Support provided by NASA Langley Research Center, through grant NAG-1-1310 is grate- 
fully appreciated. 
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Abstract— This paper investigates the buckling mode localization in the periodic multi-span beam 
with disorder occurring in an arbitrary single span. The analytical finite difference calculus is used in 
conjunction with the conventional displacement method to derive the transcendental equations from 
which buckling load is calculated. The underlying treatment is general and the solution thus obtained 
is exact. Numerical results show that the buckling mode is highly localized in the vicinity of the 
disordered span of the beam. 


1. INTRODUCTION 

Large-scale periodic structures are often encountered in engineering practice. They appear 
in different forms, such as beams on equidistant supports, gridwork structures with equally 
spaced and geometrically similar stiffening components, plates and shells with uniformly 
distributed stiffeners. Periodic structures have drawn substantial interest among research- 
ers, and many mathematical methods have been developed, starting from the pioneering 
work by Brillouin [1], The reader may also consult with studies by Lin and McDaniel [2], 
Mead [3], Sen Gupta [4], Abramovich and Elishakoff [5], Zhu et al. [6] and others. 
However, in reality, although these structures are designed to be completely identical for 
every constituent component, they seldom exhibit perfect periodicity. The deviation from 
complete periodicity is commonly known as disorder or irregularity. Disorders may arise 
from various imprecisions in the fabrication process or from geometric and material 
variations in the different parts of the structures. From the standpoint of structural analysis, 
if the constituent units of the structure are different from one another, we have to analyze 
each of them separately, with attendant satisfaction of continuity conditions from one unit 
to another. This procedure usually results in matrices of high order if the structure is 
composed of a large number of units. Owing to the fact that the inversion and other 
operations on such matrices may be involved in the calculations, the numerical errors are 
almost unavoidable. Therefore, it is often desirable to reduce the order of the matrices as 
far as possible. Fortunately, many large-scale engineering structures are essentially periodic 
and disorders often take place in some localized areas of the structure. For the analysis of 
perfectly periodic structures, the method of the finite difference calculus [7, 8] appears to 
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be very instrumental. For the treatment of periodic structures, this method has decisive 
advantages over the conventional matrix methods used in the structural analysis. The 
method usually leads to a determinative matrix that can be orders of magnitude smaller 
than those necessary in the force or displacement method, or in the conventional finite 
element method. With the size of the matrix reduced, the numerical accuracy is improved 
and the computational effort is cut down dramatically. However, the applicability of this 
method is only confined to those structures which are uniform in the spacing and stiffness 
characteristics of the constituent elements. When the periodicity is disturbed, this method 
can no longer be applied. 

Analysis of disordered structures has become of interest more recently, mainly due to the 
localization phenomenon which arises with disorders. In the past, most of the research on 
the localization phenomena have been conducted in connection with vibrations [9-12]. The 
localization phenomenon relative to buckling problem was first addressed by Pierre and 
Plaut [13], who examined the occurrence of buckling mode localization in a simple 
two-span disordered column. Ariaratnam and Xie [14] investigated the localization in the 
buckling of a system of rigid bars connected with springs, in the stochastic setting. 
Recently, Nayfeh and Hawwa [15] studied the similar problem in a more complex setting of 
multi-span columns by use of the transfer matrix method. The deterministic buckling 
localizations in cylindrical shells and in elastic plates were investigated, respectively, in 
several studies by El Naschie [16-19] and by the present authors [20]. Here, we discuss the 
general A-span beam with torsional springs at supports, which is structurally periodic 
except that one of the spans of the beam contains a disorder. By combining the finite 
difference calculus with the conventional displacement method, we present the exact 
solution for the buckling of a large-scale, multi-span periodic beam having disorder in an 
arbitrary, single span of the beam. It is shown that even a single disorder could be 
responsible for the highly localized pattern of buckling modes. 


2. BASIC EQUATIONS 


The governing differential equation for the typical ith span of the axially compressed, 
continuous beam with uniform cross-section reads 


or 



( 1 ) 

( 2 ) 


where k = V P/EI\ P is the axial load on the beam; E is the Young’s modulus, and I is 
the moment of inertia of the cross-section of the beam; y is the deflection of the beam and 
a, is the length of the ith span (Fig. 1); M?~ l , Aff are the bending moments at two 
supports of that span, respectively. The superscript l R’ (‘L’) indicates that span of the 
beam is to the right (left) of the support in question. 

The general solution to equation (2) is 


y = C t sin (fct,) + £),cos(fct i ) 4- - l) - — (3) 

P \ a, / P di 

where C, and D, are arbitrary constants which are to be determined by the use of boundary 
conditions. 
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mU 

1 — 

ka t 

Mf 

1 — 

ka t cos (ka,) 

k 2 a t EI 

J. 

sin {ka ,) . 

k 2 a t EI 

X 

sin (kai) 


We may also express Mf_ x , M~, the bending moments at supports i — 1 and i 
respectively, in terms of the rotational angles and 6 i+1 as foUows: 


where 


MU = —[icA-1 + c 2 e,] 

at 

Mi = -—[2 cA + c 2 0,_ x ] 

a, 

fot,[sin (ka t ) - ka t cos (ka t )] 
4[2 — 2cos(fca ; ) - ka, sin ( ka,)] 
ka t [ka, — sin (fca,)] 

2[2 — 2cos(^a I ) — ka t sin (ka t )] 


( 10 ) 


( 11 ) 


Note that the above derivation could also be done straightforwardly by using the stability 
function [21] . Equilibrium at a typical ith support requires that 


- M? = Jdi 


( 12 ) 


where J is the torsional modulus of the spring (Fig. 1), 6 l is the rotational angle at 
support i. 

Using equation (10), we obtain 

Ci(0i+i + 6i- i) + 4( c i + v)9 i = 0 (13) 

where v= Ja/8EI. 

Introducing the shifting operator E [8] which is defined as EQ t = 6 i+1 , equation (13) can 
be rewritten as 


[c 2 (E + E~ x ) + 4(cj + v)]0, = 0. (14) 

Equation (14) is a second-order finite difference equation, whose solution is obtained by 
letting 


6 = exp ((pi). 

A direct substitution of equation (15) into equation (14) yields 

2(ci + v) 


cosh (p = - 

Three different cases must be considered. 


C2 


(15) 

(16) 


Case A. — 2{c x + v)/c 2 5= 1. 

The solution of equation (16) has the following form 


0i,2 = ±0; 

The attendant solution for 6, reads 


0 = cosh ‘ 


2(ci + v) ' 
c 2 


6 t = Ae^ 1 + Be-* 1 

where A and B are arbitrary constants. 


(17) 

(18) 
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$1,2 = ±(ar + jtt); oc = cosh 1 
and the solution for 6 t is 


2( Cl + v) 


Cl 


j = V- 1 


0; = [A cosh (ari) + S sinh(ari)] cos (jrz). 


(19) 

( 20 ) 


Case C. — 1 «£ — 2 (c t + v)/c 2 1. 

The solution to equation (16) now becomes 


and, for 6 h we have 


$ 1,2 = ±# 


— cos 1 


2(ci + v) ' 
Cl 


Qi = A cos (/Jz) + B sin (/3i). 


( 21 ) 

( 22 ) 


Suppose that we have an A-span beam (Fig. 2), which is periodic except that its 
(q + l)th span contains a span length imperfection which makes that span slightly longer or 
shorter than the other spans of the structure, i.e. 


a, = a for i = 1, . . ., q, q + 2, . . ., N\ a 9+J = b (a A b ). (23) 


As a particular case b = a, we recover the original, perfectly periodic structure. 

To facilitate the solution of the problem, we can treat the entire beam as being 
composed of three segments. The first g-span periodic beam constitutes segment I; the 
(q + l)th span, namely, the disordered span, constitutes segment II; and the last 
(N - q — 1) spans of periodic beam represent segment III. Assume first that both segment 
I and III themselves contain a large number of spans. For segments I and III per se, the 
finite difference calculus is applicable due to their structural periodicity. As to the 
disordered span, segment II, a separate consideration should be made, and the conven- 
tional displacement method is used here. By following this procedure, we construct a 
solution composed of three parts with each part corresponding to a specific segment of the 
beam. Continuity conditions between those adjacent segments are utilized in combination 
with boundary conditions at the ends of the beam to establish an eigenvalue problem. 

For the first q spans of periodic beam, we perform the finite difference calculus analysis 


6 r = B 1 /, (r = 0, 1, . . ., q) 
where the superscript denotes the 


(24) 


cpniipncp rmmhpr nf thp 
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Fig. 2. A multi-span beam with a single disorder in the ( q + l)th span. 
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one of the three forms represented by equations (18), (20) and (22), depending on the 
physical and geometrical conditions of the segment. 

For the disordered span, recalling equation (10), we have 


ns i ” all 


M ^ +1 = - 


[2ci0o + c 2 d i 1 ] 

[2 M? + c 2 0?] 


b 

2 EI^. 


(25) 


or in another form 


where 


0 ? = 


b 2 C x M q + C 2 M q+1 


2EI 


4c i - cl 


0? = - 


b 2diM q+ i + c 2 M q 


2EI 


4c? - c 2 2 


Cl 


fcb[sin (kb) — kb cos (kb)] 
4[2 - 2 cos (kb) — kb sin (kb)] 
kb[kb — sin (kb)] 


(26) 


(27) 


2[2 - 2 cos (kb) — kb sin (kb)] 

Note that c 1 and c 2 are obtained from equation (11) by formally replacing a t with b. 

The treatment of the last N - q - 1 spans of periodic beam is similar to that of 
segment I, 


0 , = 0 


m 


-q- 1 ’ 


(s = q + 1, q + 2, . . ., N) 


(28) 


where 0;_ 9 _i again adopts one of the three forms denoted by equations (18), (20) and (22). 

Consider now a beam simply supported at its two ends (other boundary conditions can 
be treated in a similar manner). Then the boundary condition at the left end of the beam 
can be represented as 

-Mg = JO o or (2c x + 4 v)0q + c 2 0\ = 0 (29) 

while the boundary condition at the right end of the beam reads 


Mi 


T() m 

J VN-q - 1 


or 


(2c I + 4 v)0™_ 9 _! + C 2 0^_ 9 _ 2 


3 in 


0. 


(30) 


Conditions of continuity between the periodic spans and the disordered span of the 
beam, namely, between segment I and segment II, are 


Mi 


Mq = JOq 


0'q = 0 ? 


or (2c 1 + 4 v)0q + c 2 0i 1 + 


-M R q 


0 


or 




— Z—Aui — 

a(4c\ — cl) \ 2 E 


2EI 


Mq + C 2 — — 

2 El 2EI 


= 0. 


(31) 


Analogously, the continuity conditions between the second and the third segments are 

or (2c ! + 4v)0o" + c,df - —ML , = 0 
b 


M L q+l - M q+1 ■■ 


J0 l o" 


2EI 


0 ? = 0 " 


or 


0o" + 


a(4c? 


cl) 




M* + 2c 1 


11 

Mq+l 


2 El 




( 32 ) 
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Equations (31) and (32) should be formulated in terms of the three different cases, 
because in each case, the resulting expressions are different. 

For Case A, the rotation angles in the first and third segments can be expressed as 

6\ = A&+ + B&-+ (r = 0, 1, . . ., q) 


0 


III 

s—q — 1 


AiQ M.s- q -\) + Bi& -<Hs-q- 1) (s = q + 1, . . N). 


Substituting the above expressions in the boundary conditions (29) and (30) and the 
continuity conditions (31) and (32), we obtain six homogeneous algebraic equations, 

(2cj + 4v + c 2 e 4 ’)Ai + (2c 1 + 4v + c 2 e~‘ t> ) B 1 = 0 


[(2 c x + 4v)e*« + cze^- 1 ^! + [(2c x + 4v)e“* + + M* = 0 


(34) 

(35) 


e H Al + e ~* Bl - 2ci J -)m« - t C2 - (-)m^ = 0 

4c\ - c\\a ) 4cf - c$\a J 
(2 Cj + 4v + c 2 e 4 ’)A 2 + (2 Cj + 4v + c 2 e~' p )B 2 - M^ +1 = 0 


A 2 + B 9 + 


Cl 


4 cj + cl’ 


( b\o7 R i ^Ci / t>\xi L _ 
I M, + - 2 I l-^g+1 

V a ] 4cj — c 2 \a ) 


0 


[(2cj + 4v)e' KN " ? - 1) + c 2 e^ N - q ~ 2) ]A 2 + [(2c x + 4v)&~^ N ^ l) + c 2 &-^ N - q ~ 2) ]B : 
where 

- M " a - »„ - 




fL _ Mq + i a 


2 El 2EI 

For Case B, the solutions for the first and third segments are as below. 

Q\ = A x cosh (ar) cos (irr) + B x sinh (ar) cos (nr) (r = 0, 1, . . ., q ) 


(36) 

(37) 

(38) 

2-0 (39) 
(40) 


6 


hi 


(41) 


(s = q + 1, 


N) 


7s-q - 1 A 2 cosh[ar(s - q - 1 )]cos[jt(s - q - 1)] 

+ B 2 sinhfofs — q - 1)] cos [77(5 - q - 1)] 

and performing the substitution similar to that in Case A, we arrive at the following six 
homogeneous equations 

[2 Cj + 4v - c 2 cosh (a)]A l - c 2 sinh (a)B 1 = 0 (42) 

{(2c! + 4v) cosh (aq) cos (nq) + c 2 cosh[ar(g - 1)] cos [tt (<7 - 1)]}A X 

{(2 Ci + 4v) sinh (aq) cos (nq) + c 2 sinh [a(q - 1)] cos [n(q — !)]}/?[ + M q = 0 (43) 


2cj ( b \ — r c 2 ( b\ 

cosh (aq) cos (nq)A 2 + sinh (aq) cos (nq)B x -I — \M q - — -I — 

4c\- c\\ a) 4c{- c\\ a } 


Mq+l = 0 


[2 Ci + 4v — c 2 cosh(ar)]A 2 — c 2 sinh(a)B 2 - M q+l = 0 

a, + + ~~i (-) - 0 

4 c\- c\\ a) 4 c\- c 2 \ a] 

{(2ci + 4v) cosh [a(N - q - 1)] cos [7r ( N - q - 1)] 

+ c 2 cosh[ar(iV - q - 2)] cos [n(N — q — 2)]}A 2 
X {(2cj + 4v)sinh[ar(A - q - 1 )]cos[ 7 t(A — q — 1)] 
+ c 2 sinh [cr(iV — q — 2)] cos [tt(N — q - 2)]} B 2 = 0. 


(44) 

(45) 

(46) 


( 47 ) 
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For Case C, the solutions are in the following form 
=A 2 cos(pr) + fijsin()3r) 

= A 2 cos[P(s - q - 1)] + B 2 sin [/?(s - q - 1)] 
and the corresponding equations are 

[2 Cx + 4v + c 2 cos(/3)]y4 1 + c 2 sm(P)B t 
{(2c x + 4v)cos (Pq) + c 2 cos[P(q - 1)]}Aj 


(r = 0, 1 , . . q) 

(s = q + 1, . . N) 


(48) 


(49) 


+ {(2cj + 4v)sin(/?g) + c 2 sin[P{q — l)]}Si + Mq = 0 (50) 


cos (Pq)A 1 + sin (pq)B 1 j (— W* - 2 J — )m^ +1 

4c\ - c 2 \ a ) 4c\ - c 2 \a / 

[2ci + 4v + c 2 cos (P)] A 2 + c 2 sm(P)B 2 - M^ +1 = 0 


0 


A-y + 


+ -#-r (-) 

4cj - c 2 \a) 4c\ - c\ \ a ) 


Mq+\ = 0 


4ci c 2 

{(2c 1 + 4v)cos[P(N — q — 1)] + c 2 cos[/?(iV — q — 2)]}A 2 

+ {(2ci + 4v)sin[/?(iV - q — 1)] + c 2 sin[P(N — q — 2)]}S 2 = 0. 


(51) 

(52) 

(53) 

(54) 


Thus, for each case, we have six homogeneous algebraic equations, which can be 
expressed in a matrix form as follows 

[F(K)] 6x 6 {5} 6x1 = 0 (55) 

where [F(f£)] is the coefficient matrix, K = ka and {b} 3 * * * 7 = {A u B r , Mq, Mq +1 , A 2 , B 2 }. 

Non-triviality of {6} requires that the determinant of the coefficient matrix vanish, 

Det[F(/C)] = 0 (56) 

which constitutes a transcendental equation for the non-dimensional buckling load para- 
meter K. Once the buckling load parameter K is determined, we can use equation (6) to 
calculate, span by span, the buckling mode shapes for the entire structure. 

It is worth mentioning that, if the disorder occurs in the first or last span of the beam, 
the whole beam can be partitioned into two segments; one is the disordered span and the 
other is the N- 1 spans of periodic beam. If this happens, only four equations are needed 
to characterize the problem so that instead of having a 6 x 6 matrix for [F(/sT)], we will 
have a 4 x 4 deter m inant matrix (see Appendix for details). 


3. EXAMPLES AND DISCUSSION 

In the following numerical examples, the non-dimensional spring constant v is fixed at 
0.3 since this particular case for perfectly periodic beam was considered by Wah and 
Calcote [8] and a comparison can be made with their results. 

As a first example, we discuss a simply supported, 100-span continuous beam. The 

disorder arises from a span length inperfection characterized by b/a = 1.1, where a and b 

are the lengths of periodic spans and the disordered span, respectively. Since b/a > 1, the 
disordered span is longer than other spans of the structure. The disorder may appear in any 
span of the beam. Figure 3 depicts the results of the buckling load parameter K for the 
beams with or without torsional springs. The most critical situation for the beam without 
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Fig. 3. Change of buckling load with the location of the disorder (i— the sequence number of the span where the 

disorder occurs). 


torsional springs occurs when the disorder occurs at either the first or the last span, for 
which the buckling load parameter K equals 3.05. If the disorder appears in one of the 
spans close to the center of the beam, the buckling load increases. However, for beams 
with torsional springs, numerical results display a quite different picture. For this case, the 
occurrence of disorder near the boundaries may be more advantageous, since the buckling 
load decreases as the disorder moves away from the boundaries. Nevertheless, for both 
cases, the buckling load remains almost unchanged with the location of disorder, once the 
disorder is 10 spans away from both boundaries. This means that the effect of the boundary 
dies out if the disorder is sufficiently far from the boundary. The location of the disorder 
has almost no effect on the buckling load if the spring modulus v = 0.15. If we refer to the 
buckling load in the absence of disorder as the classical buckling load, then the classical 
buckling load parameter K equals it for the case without torsional springs, and has a value 
of 3.760 for the beam with torsional springs of v = 0.3 [8]. With the disorder present, 
numerical results show that the buckling load parameter K is below the corresponding 
classical value. For instance, K is less than 3.12 for the beam without torsional springs, and 
is no more than 3.72 for the beam with torsional springs of v= 0.3. Thus, we can see that 
such a disorder, namely, the span length imperfection, may have a degrading effect on the 
load-carrying capacity of the structure. 

The second example is mainly devoted to the discussion of the buckling mode shapes for 
disordered periodic beam. Figures 4 and 5 depict the buckling mode shapes for an 11-span 
beam with the disorder appearing at different locations of the beam. Again, the disorder is 
introduced by a span length imperfection specified by b/a — 1.1. A significant phenomenon 
is that the buckling mode shapes exhibit a strong localization in the disordered span, when 
the torsional springs are used at supports. The larger the moduli of the torsional springs, 
the more localized the mode shape becomes. Thus, the torsional spring weakens the 
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Fig. 4. Buckling mode shapes for a disordered 11-span beam without torsional spring (i— support sequence 
number, y— deflection), (a) Disorder occurs in the third span; (b) disorder occurs in the mid-span (sixth span); 

(c) disorder occurs in the last span. 


coupling between different spans of the structure. This observation is consistent with that 
found by Pierre and Plaut [13] for the two-span beam. 

In passing, it is worthwhile to point out that, although the underlying treatment makes it 
possible to obtain an exact solution to the buckling problem of disordered periodic beams 
with any number of spans by dealing with a determinative matrix of very low order (the 
matrix is 6 x 6 if only one disorder occurs in the span neither the first nor the last), some 
numerical problem can occur when N, the number of spans, is a large number, say 
IV 3= 50. This is because, in our calculations, we have to evaluate terms e^ N ~ q ~ 1 \ 
cosh [<v(/V — q — 1)], which can be so large when q is small, that they may exceed the 
upper limit of a digital computer. To avoid this numerical difficulty, we may divide 
the corresponding equation by a relevant, large term, for example, e^N-g-i) or 
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Fig. 5. Buckling mode shapes for a disordered 11-span beam with torsional spring of v= 10 (i— support sequence 
number, y— deflection), (a) Disorder occurs in the third span; (b) disorder occurs in the mid-span (sixth span); 

(c) disorder occurs in the last span. 


cosh[<*((V - q - 1)], and manipulate the resulting equation by making an asymptotic 
approximation 

sinh [ft(JV - q - 1)] _ , , « N , ( 57 ) 

cosh[<*(./V — q — 1 )] 

It follows that only equations (39) and (47) need to be modified and they adopt, after the 
approximations, the following forms, respectively: 

[(2cj + 4v) + c 2 e-^]A 2 = 0 (58) 

and 


{(2c x + 4v)cos[tt(N — q — \)\ + c 2 e a cosjy(iV — q — 2)]}A 2 

+ {(2c x 4- 4v)cos [tt(N -<?-!)] + c 2 e _ “ cos [n{ N — q — 2 )]}B 2 = 0. (59) 
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Figure 6 depicts the buckling modes of a 100-span beam and a 400-span beam with 
torsional springs of v= 0.3; both beams contain a disorder in the 40th span. From Fig. 6, it 
is clear that the deflection of the beam at buckhng dies out quickly as the distance from the 
disordered span increases. The envelope of the buckling mode is depicted in Fig. 7. If we 
take a logarithm of the function, we obtain Fig. 8, which displays a nearly straight line. 
Thus we establish that the deflection at buckhng decays exponentially. The exponential 
decay constant [12], sometimes referred to as the Lyapunov exponent [14] in the literature, 
(its counterpart in vibration problems is commonly known as the logarithmic decrement 
[22]) equals -0.260. 

In the above examples, we have dealt with the disorder of ‘detrimental’ character; that 
is, the disordered span is longer than the periodic spans, and such a disorder leads to a 
reduction in the buckling load of the structure. As a third example, we consider also 
another possibility with attendant opposite effect; that is, the disordered span being shorter 
than the other spans ( b/a < 1). In this case, the disorder turns out to be ‘beneficial’ 
because it results in slightly higher buckling loads. For example, for a perfectly periodic 
11-span beam with torsional springs of non-dimensional modulus v= 0.3, the buckhng load 
parameter K is, as mentioned earlier, 3.76; for the same structure but with disorder, K 
varies in the range from 3.79 to 3.82, depending on the location of the disorder. Thus, we 




Fig. 6. Buckling mode localization in multi-span beams with torsional springs of v= 0.3 (the disorder occurs in the 
40th span), (a) Buckling mode for a 100-span beam; (b) buckhng mode of a 400-span beam. 
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Fig. 7. Envelope of the buckling mode shape for a 100-span beam. 



Fig. 8. Logarithmic plot of the envelope function. 



Fig. 9.(a). 
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Fig. 9. Buckling mode shapes for a disorder 11-span beam (disordered span is shorter than the periodic span), 
(a) Disorder occurs in the third span; (b) disorder occurs in the mid-span (sixth span). 


can see that the ‘beneficial’ effect of such a disorder on the buckling load is very small— the 
increase in buckling load is less than 2% . However, the impact of disorder on the buckling 
mode shape is noteworthy. Figure 9 portrays the buckling mode shapes of 11 -span beams 
containing such a disorder, from which it can be seen that the change in buckling mode is 
significant, especially when the disorder occurs near the boundaries. Interestingly enough, 
the localization phenomenon has not been obseved in this ‘beneficial’ case. 

Finally, it is should be noted that the present study can be generalized to other more 
complicated problems. The work is underway and will be reported elsewhere. 
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APPENDIX 


Consider the case when the disorder occurs in the first span. We have the following set of equations for 
boundary and continuity conditions. 

(1) Boundary condition at the left end 


r—T + 1 M 0 R + J-) Cl Mi = 0. 
i - cl \ a ] 4cf - c\ 

(Al) 

Mi - (2 Cl + 4 v)8y - c 2 8{ ! = 0. 

(A2) 

-(c 2 M 0 r + 2 CiMt) = 0. 
a(4c\ - ci) 

(A3) 

(2d + 4v)0#_ 1 + c 2 6n- 2 = 0. 

(A4) 


-M$ = J6o or 

(2) Equilibrium equation for support 1 

Mi -Mi = jetf or 

(3) Continuity condition between segments I and II 

e{ = ey or ey + 

(4) Boundary condition at the right end 


Again, three different cases have to be considered separately. The expressions for 0 ! J in those three cases are as 
follows. 


Case A: 6 1 / — A 2 c' t ‘ s + B 2 e~' l ’ s 

Case B: 8jf = A 2 cosh (ars) cos (?is) + B 2 sinh (as) cos (jts) (A5) 

Case C: 8'/ = A 2 cos(/5j) + B 2 sin (fis) (s = 0, 1, 2, . . ., N — 1). 

Substituting equation (A5) into equations (A1)-(A4), one obtains four homogeneous algebraic equations for 
each case. These four equations can be expressed in a matrix form 

[F(K)] ix4 {6} M = 0, (A6) 


where [F(A)] is a coefficient matrix, and {8} T = {M§, Mf , A 2 , B 2 ). 

It is of interest to note that the above equations can be formally obtained from equation (55) by letting q = 0 
and Ai = Bi- 0. As to the case where the disorder occurs in the last span of the beam, a similar treatment can be 
made accordingly and the problem again reduces to a set of four homogeneous algebraic equations. 
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1 Introduction 

Composite materials are widely used in various types of engineering structures. 
Previous studies have been based on the assumption that the properties of the 
composite materials are characterized by pre-determined elastic moduli, and no 
uncertainties of these moduli have been considered (Vinson and Sierakowski [1986], 
Tennyson [1975], Whitney [1987]). However, the composite materials are subject 
to a certain amount of scatter in their measured elastic moduli (Tewary [1978]). 
Such uncertainties in elastic moduli are due to many factors which influence the 
actual properties of composites. For example, among other things, misalignment 
of fibers or imperfect bonding between the fibers and the matrix may contribute 
to the scattered values of the measured elastic moduli. To a large extent, the 
properties of composite materials are dependent on the fabrication process. But 
even the composite materials manufactured by the same process may demonstrate 
differences in their elastic properties. For design purposes, one should be aware 
of the potential variations in load-carrying capacity and dynamic behavior of such 
structures that can arise due to the uncertainty in elastic moduli. A more realistic 
analysis of composite structures should be performed with the variations of the 
elastic moduli being taken into consideration at the same time. 

Scatter in material properties is usually treated within the realm of the stochas- 
tic finite element method (Anantha and Ganesan [1993]). However, as Shinozuka 
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[1987] mentions: “...it is recognized that it is rather difficult to estimate experimen- 
tally the auto-correlation function, or in the case of weak homogeneity, the spectral 
density function of the stochastic variation of material properties. In view of this, 
the upper bound results are particularly important, since the bounds derived ... do 
not require knowledge of the autocorrelation function.” In the stochastic context, 
Shinozuka [1987], and Shinozuka and Deodatis [1988] and Deodatis and Shinozuka 
[1989] derived upper and lower bounds on the probabilistic characteristics of the 
response in terms of probabilistic characteristics of the material variability. 

Recently, in a private correspondence, Ariaratnam [January, 1992] advocated 
placing bounds on elastic moduli and deriving appropriate bounds for certain struc- 
tural properties such as the buckling loads and natural frequencies. This idea is 
implemented in this investigation. Here, we use convex modeling (Ben-Haim and El- 
ishakoff [1990]), a newly-developed analytical tool, which incorporates uncertainties 
in elastic moduli into the structural analysis. Since it came into being, convex mod- 
eling has been used under different circumstances for solving a variety of engineering 
problems (Elishakoff and Ben-Haim [1990], Elishakoff, Li and Starnes [1994]). 

The present paper is a generalization of the previous study of Elishakoff, Li 
and Starnes [1994], where the influence of uncertainty in elastic moduli on the 
axial buckling load was discussed. Here, we consider another case of buckling, that 
is, shells under uniform external pressure. In addition, this paper deals with the 
variability of natural frequencies by use of convex modeling, which is apparently the 
first study of this kind in the literature. A numerical approach to the uncertainty 
problem is nonlinear programming, which we apply to solve the same problem 
to generate a set of comparable numerical data. The results from both methods 
show good agreement throughout. Thus, the effectiveness of the analytic convex 
modeling is clearly demonstrated. The bounds of the natural frequency and the 
buckling load provide the designer with a better view of the vibrational behavior 
and the actual load carrying capacities possessed by the composite structure. 


2 Basic Equations 

We use the Donnell shell theory (Timoshenko and Gere [1961]) for the analysis 
of buckling of cylindrical shells of composite materials. The strain-displacement 
relations are 


€-X 

du 
dx ’ 

f^x — 

6y = 

dv | w 
dy IV 

Ky = 

7xy = 

dv i du 
dx ' dy ’ 

foxy ~ 


d 2 w 
dx 2 
d 2 w 
dy 2 
2 d 2 w 
dxdy 


( 2 . 1 ) 


where x and y are the axial and circumferential coordinates in the shell middle 
surface; u and v are the shell displacement along axial and circumferential direc- 
tions, and w is the radial displacement, positive outward; e x , e y and j xy are strain 
components; k x , n y and n xy are middle surface curvatures of the shell; R is the 
radius of the cylindrical shell. The constitutive relations for the composite laminate 
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read 
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^4i6 A26 Aqq B i6 B26 Bee 


e xy 
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B 11 B 12 i ?26 T>h D 12 Die 


K>x 
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B12 B22 B26 D 12 D22 £26 


Ky 

\MxyJ 


\-Bi6 B26 Bee Die T>26 De6 / 


\K xy J 


where N x , N y and N xy are stress resultants, M x , M y and M xy are bending and 
twisting moments, acting on a laminate; the laminate stiffnesses A- t j, Bij and Dij 
are defined as 


fh/2 

B^, Dij) = / Q\ 3 \\,z,z 2 )d. 

J-h/2 


where h is the total thickness of the laminate, and 2 is the coordinate in the di- 
rection of the laminate thickness; Qy are the transformed reduced stiffnesses and 
can be expressed in terms of the lamina orientation and four independent engineer- 
ing material constants in principal material directions, i.e., E\, E 2 , ^12 and G 12 
Jones [1975]. 

The equilibrium equations of the cylindrical shell read 


dN x dN xy _ 
dx dy 
9N xy dN y 
dx dy 

d 2 M x d 2 M xy d 2 M y 

dx 2 dxdy dy 2 


(2.3) 


R Ny + Nx dx 2 +2Nxy dxdy 


d 2 w d 2 w 
+ N V 


dy 2 


-p 


d 2 w 

W 


= 0 


where p is the mass per unit volume of the shell and t is the time variable. 
Using equations (2.1) and (2.2), equation (2.3) can be written as 


( T 11 T 12 T i3 \ / u\ 

L21 L22 T23 I v ) 
T31 T32 T33/ \wj 


0 

0 

d 2 w 


Q^2 + 2N xy dx Q y + N r 


d i W 

y dy 2 



(2.4) 
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where the operators L i3 are 


Lu = 

L\2 = 

-hl3 = 


L 22 = 
L 23 — 


£33 = 


. d n . d 2 . d 2 
A U + 2Ai 6 • + Aqq—-:- 

dx 2 dxdy ay 2 

d 2 () 2 {)- 

Aie dx* + (Al2 + Am) forty + A ™W 

1 , d_ d_ d 2 

R^ Al 2 dx + 26 dy^ Bn dx 2 


— 3Si 6 


d 3 


dx 2 dy 


d 3 d 3 

{Bl2 + B66) fort^~ B26 W 


. d 2 nA d 2 

A66 fo^ +2A26 forty 
R^ 26 dx +A ‘ 22 dy ) ' 


. d 2 
+ A22 df 

dy 3 


B 22 


-m 


d 3 


26 : 


dxdy 2 
- — + 2B 2 6 


(B\2 + 2Bqq) 
d 2 


d 3 


R 
+4Di6 


' dx 2 

d 4 

dx 3 dy 


dxdy 


■ 2{D\2 + 2 £) 66 ) 


dx 2 dy 
d 2 

+ B22 W ) + 

d 4 

dx 2 dy 2 


B 16 


dx 3 

A 22, d 94 

~W +Dll w 

d 4 


+ 4£>26 


(2.5) 


0 4 


dxdy 3 ”^ 22 dy A 


We consider the cylindrical shell with simply supported boundary conditions 
which are satisfied by the following displacement functions, 


w\ co 00 /UmnCos(X m x)cos(X n y)e tuJt \ 
v ) = 11, I Vmn sin ( A mtf) sm(X n y)e iu,t I (2.6) 

wj m=i«=o \W mn sm(X m x)cos(X n y)e iait ) 


where A m = rrnr/L, X n = n/R. and u> is the natural frequency of the shell. 

Similar to Tasi [1966] and Hirano [1979], here the coupling stiffnesses 
(Aie, A 2 q, Bie, B 2 6, Die, D 26 ) are neglected. They actually vanish for symmetric 
cross-ply laminates. As for symmetric angle-ply laminates, B\e and B 2 6 are zero, 
and Am, A 2 6, T>i6 and D 2 e can be neglected for laminates with many layers. 

Substitution of equations (2.1) and (2.2) into equation (2.4) leads to a set of 
homogeneous linear algebraic equations, and the existence of non-trivial solutions 
requires that the determinant of the coefficient matrix vanish, 

fC n C12 C13 \ 

det | C 21 C 22 C 23 ) = 0 (2.7) 

\C .31 C 32 C 33 — A“ n N x — X rn X n N xy — X^N y + pin 2 } 
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where elements Cij ’s are expressed as 
Cn = 

C22 = A 22 X^ + Aq 6 A^j 

C33 = -DnA^j + 2 (Di 2 + 2 £> 66 )A^A^ + £> 22 A^ + + 2 -^-A^ + 2 ^^-A^ 

C12 = C21 = (-A12 + ^ 66 )^mA re ( 2 . 8 ) 

C23 = C32 = (B12 + 2 i? 66 )A^ t A rl + — ^~A n + B 22 X'^ l 

C13 = C 13 = ^A m + BnA^ + (B 12 + 2 B e 6 )X m Xl 

IX 


Prom equation (2.7), the following expression can be readily derived for the 
natural frequency 


puj 


2 — C33+ 


(2C12C23 - C' 13 C 2 2)C 1 3 - C&Cu 


C ?2 


■ C11C22 


(2.9) 


+ A m A(r + A m A n N xy + X n N y — puj. 


mn,N 


where the subscript “N” indicates the presence of external loading acting in the 
mid-surface of the shell. 

If the shell in question is free from external loading, i.e. N x = N y = N xy = 0, 
the natural frequency becomes 


1L . (2C12C23 - C 13 C 22 )C n - C| 3 Cn 
p L° 33 + (C? 2 - C n C 22 ) 


(2.10) 


where the subscript zero indicates that the external loads acting in the mid-surface 
of the shell are absent. 

To determine the fundamental natural frequency for a cylindrical shell with 
given dimensions and material properties, one determines those integer values of m 
and n which minimize uj mn . 

If w = 0, equation (2.9) yields an expression for the buckling load. The case of 
static axial buckling of composite shells has been discussed in our previous study 
Elishakoff, Li and Starnes [1994] . Here we will confine ourselves to the investigation 
of the buckling of shells under external pressure p, for which 


Ny = -pR, N x 


pR 

T’ 


N X y =0 


( 2 . 11 ) 


The expression for the critical external pressure can be readily derived as follows 


P = 



C33 + 


(2C 12 C 23 - C i:i C 22 )C 13 - CfgCn 
Ci 2 2 - C11C22 


= Pmn (2.12) 


Again, one has to perform a search with respect to integer variables m and n to 
minimize the objective function p mn in order to obtain the critical external pressure 
p cr for a cylindrical shell with given dimensions and material properties. 
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3 Convex Modeling of Uncertain Moduli 

As has been stated in the previous section, the objective function (the natural 
frequency w m „ i0 or the critical external pressure p mn ) of the structure depends on 
the four basic elastic moduli. (One can also discuss the natural frequency uj m n,N 
if the external pressure is present, provided that the external pressure is below 
its critical value.) For the sake of generality, in the following analysis, a generic 
formula for the objective function is adopted instead of relying on a more concrete 
expression such as equation (2.10) or (2.12). 

The objective function is written in the following generic form, 


£ — £(£ l ? £2, V 12, U12), £ — Wmn, 0 Or Pm n (3.1) 

or more simply, 


. F = F(Ei), (i = 1,2, 3,4) (3.2) 

where £3 = iq 2 and £4 = Gi2- The function £ in the above equation also depends 
on the form of structure, boundary conditions as well as geometric properties. 

Let Ef(i = 1, 2, 3, 4) be the nominal values of the elastic moduli, which might be 
visualized as the nominal values of the elastic moduli, or the average values derived 
from test data. Then, the elastic moduli of values different from those nominal 
values could be denoted as Ef + S if Si being deviations from nominal values. The 
objective function corresponding to these elastic moduli, retaining only the first 
order terms in Si, is written as follows: 

F(Ef + Si) — F(Ef) + f T S (3.3) 


where 

,T ,, , , , , ,T (0FIF»i SF(E ° ) SF(Ef) dF(E», \ 

5 -(*!.*>.«».«. / -{-dET'-dEr’^gEr-SEr) ( ’ 

The deviation S from the nominal elastic moduli is assumed to vary in the 
following ellipsoidal set: 


Z(a,e) = {S: £ f < a 2 } (3.5) 

i= 1 e * 

where the size parameter a and the lengths of semi-axes e, (i = 1,2,3, 4) of the 
ellipsoid are based on the experimental data available, and will be discussed later. 

The problem is formulated as follows: given an ellipsoid of the elastic mod- 
uli (equation (3.5)), find the extremal natural frequency (or the critical external 
pressure) 

F ext = extremum 
6eZ( a ,e)(F(E°) + f T &). 

In equation (3.6), F ext is the lowest or the highest value of the fundamental 
natural frequency (or the critical external pressure) of the composite structure with 
the elastic moduli varying within the range of the ellipsoidal set Z. Since equation 
(3.6) calls for finding the extremum of the linear functional f T 6 on the convex set 
Z(a, e), the extremal values take place on the set of the extreme points, or the 
boundary, of the set Z as discussed by Ben-Haim and Elishakoff [1990] . 
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To arrive at the extreme values, we use the method of Lagrange multipliers. 
Since the analysis which follows is mathematically analogous to that described in 
the monograph Ben-Haim and Elishakoff [1990], here we only list the final result. 
The extremal values of the objective function F result in the following expressions 


LJ. 


mn, max 


CO. 


mn, min 


— < ^mn, 0 ^ 


— ^ ^mn,0 ^ 


1 


Cj duJjnn^o^Ef) ^2 


^ < ^mn,0 > 9E, 


1 - a. 


V~W e i 9ui m n,o(Ef) ,2 

Aj < u mn> o > dEi 


Per , max — ^ Per ^ 


Per, min 


^ Per ^ 


1 + O' A 


1 — a 


4 

V'V e * 

9Pcr(Ei)}2 


dEi > \ 

4 

dPcr(Ef). 2 

^<Pcr> 

dEi > \ 


(3.7) 


where <C ^mn,o ^ — ■u’mn.o (F^i , T']) , Eg , T) 1 ) and < p cr > — 7A:r(L'i , 7^2 , Eg - L' l ) are 
values of the natural frequency and critical external pressure calculated at the 
average values of the elastic moduli. 

From the above equation, the upper and lower bounds of the natural frequency 
ui and the critical external pressure p cr can be calculated if we use the appropriate 
expression for F. Equation (3.7) shows explicitly that the uncertainties in elastic 
moduli have a direct effect on the values of the natural frequency and the buckling 
load of the structure. It is seen here that the products of lengths of semi-axes of 
the uncertainty ellipsoid and the sensitivity derivatives play an important role in 
the variation profile of the objective function in the convex modeling. 


4 Determination of Convex Set from Measured Data 

Before any prediction can be made on the vibratory properties of the composite 
structure, the values of a and e»(i = 1, 2, 3, 4) should be known in advance. In fact, 
these values are dependent on the manufacturing process by which the composite 
structure has been fabricated. It is understandable that the more advanced the 
manufacturing process and the better the workmanship, the smaller these param- 
eters will be in value. Besides, the evaluation of these parameters is also linked 
with the amount of information one has about the properties of the composite 
considered. In this study, the parameter a is set equal to unity. As long as a is 
fixed, the other parameters could be readily determined by the evaluation of the 
existing experimental data. Normally, if a sufficient amount of experimental data 
is available, the average value of these data could be used as the nominal value Ef 
for the corresponding elastic modulus, whereas parameters e^’s could be chosen as 
the proper deviations from the average values of the corresponding measured data. 
Even if one has only limited knowledge of the uncertain elastic moduli, it is still 
possible to make judgement on what value the natural frequency of the composite 
structure will be, provided the variation ranges of these elastic moduli are known. 
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Suppose that from the experimental data we deduce that the elastic moduli are 
varying in the following ranges: 

Et<Ei<E f (i = 1,2, 3, 4) (4.1) 

where E 1 / and Ef correspond to upper and lower bounds of the elastic modules E t 
respectively. By using the method described in Elishakoff, Li and Starnes [1994], 
the nominal values of the elastic moduli and the semi-axes of the ellipsoid can be 
determined respectively as 

&i={E¥ +E?)/2, e i = Ej f -Ef', (i = 1,2, 3, 4) (4.2) 

Thus, once the experimental data are available, the semi-axes of the ellipsoid 
can be determined, and the analysis of convex modeling can be carried out. 

5 Numerical Analysis by Nonlinear Programming 

The above problem can also be treated as a nonlinear programming problem 
with bounds, which is mathematically stated as follows: 

Find m.inF(Ei,E 2 ,E 3 ,E 4 ) or maxF(Ei,E 2 ,E 3 ,E 4 ,) (5.1) 

subject to Et 1 < Ei < E ^ , for i = 1,2, 3, 4 


For this problem, it is not advisable to apply the gradient methods or Davidon- 
Fletcher-Powell method (Himmelblau [1972]) which are used most often in nonlin- 
ear optimization, since the directional derivative of the objective function F can 
not be easily calculated analytically. So, a direct search should be implemented. 
Here we choose the complex method (Himmelblau [1972] , Beveridge and Schechter 
[1970]) which is based exclusively on function comparison and no derivatives are 
used. In the search for a minimun F, the complex method starts with 2n (n = 4 in 
our case) points E^, . . . , E^ 2 "), where EW = {E^\e^\ E^\e{^}. At each 

search cycle, a new point is generated by a certain rule in terms of the previous 
2 n points and the worst point E^-b , which has the largest value of F among these 
2n points, is rejected and replaced by the new point. Whenever the new point 
generated is beyond the bound, it will be set to the bound. Progress will continue 
with repeated rejection and regeneration until some criteria is met. For a complete 
description of this method, one may consult Himmelblau [1972], Beveridge and 
Schechter [1970]. Nowadays, with many commercial software packages available 
for numerical analysis, such as IMSL [1989], performing non-linear programming 
could also be realized through use of such computational tools as gradient projec- 
tion, feasible direction and penalty-function methods (Kirsh [1981]). As was shown 
above, since the expression of the natural frequency and buckling load is available 
explicitly, one can choose to optimize a two-term or three-term Taylor expression 
(Ben-Haim and Elishakoff [1990]), of the natural frequency and the buckling load 
about the nominal values of elastic moduli subject to nonlinear quadratic constraint 
function. The constraint function, given in equation (3.5), represents the equation 
of the ellipsoid of minimum volume that encloses the rectangular parallelepiped 
representing the original inequality constraints. 

It appears remarkable that there is a basic philosophical difference with the 
classical optimization studies, whereas in classical optimum design of structures one 
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looks for the maximization of the buckling loads and natural frequencies, here we 
look for the least favorable scenarios; namely, we determine, for the design purposes, 
minimum natural frequencies and minimum buckling loads. This procedure has 
been dubbed in Elishakoff [1991a] “anti- optimization”. 


6 Numerical Example and Discussion 


We now proceed to investigate several cases of the free vibration and buckling 
problem of composite shells with a view of gaining some insight into the effect 
of uncertainties in the material properties on the load-carrying capacity and the 
dynamic behavior of these structures. 

In the following analysis, elastic moduli used are from real material tests 
(Goggin [1973]). The material of the lamina is composed of carbon fibers and 
resin matrix, with a volume fraction of fiber being 40%. From the figures in 
Goggin [1973] , the following data were deduced 


= 100GPa(14.5 x 10 6 psi), 
E$ = 7.3GPa(1.06 x 10 6 psi), 
= 0.28, 

E% = 3.2GPa(.46 x 10 6 psi), 


E f = 90GPa(13.0 x 10 6 psi) 
E% = 6.8GPa(.99 x 10 6 psi) 
E% = 0.22 

Ef = 2.6GPa(0.38 x 10 6 psi) 


( 6 . 1 ) 


From these data, the nominal elastic moduli Ef and the semi-axes a can be 
evaluated, by use of equation (4.2), as the following 


£? = 13.75 x 10 6 psi, 
E5 = 1.03 x 10 6 psi, 
E§ = 0.25, 

El = 0.42 x 10 6 psi, 


ei = 1.5 x 10 6 psi 
e 2 = 0.07 x 10 6 psi 
— 0.06 

e± = 0.08 x 10 6 psi 


( 6 . 2 ) 


The natural frequency and the critical external pressure of the composite cylin- 
drical shell are given be equation (2.10) and equation (2.12), respectively. It should 
be noted here that the basic elastic moduli Ei are implicitly contained in the flexural 
stiffnesses Aij and Dij. 

The shells investigated have a radius 6.0 inch, are composed of 0.012-inch thick 
layers. The following two cases are considered: 

Case 1: the 10-layer laminated shell, with ply angle being [9, —9, 9 , —9, 9, — 9) sym , 
9 ranging from 0° to 90° . 

Case 2: the 5-layer laminated shell, with ply angle being [9, — 9, 9i/2] sym , 9 
ranging from 0° to 90° . 

In order to quantify the variability of the natural frequency or the critical 
external pressure of the composite structure, a percentage variability /3 is defined 
as follows 

f3 = F \~ Fl x 100% (6.3) 

2.r n 

where subscripts u, l and n denote the upper-bound, lower-bound, and the nominal 
value, respectively. 

Figures 1 and 2 portray the variability of . the fundamental natural frequency 
due to the uncertainty in elastic moduli for the above two cases. Figures 3 and 
4 depict the variation profile for the critical external pressure of a set of 10-layer 
and 5-layer laminated shells. Note that the abrupt turns at some points of these 
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Figure 1 Variability of fundamental natural frequency for a set of 10-layer 
laminated cylindrical shells 



Figure 2 Variability of fundamental natural frequency for a set of 5-layer 
laminated cylindrical shells 
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Figure 3 Variability of critical external pressure for a set of 10-layer laminated 
cylindrical shells 



Figure 4 Variability of critical external pressure for a set of 5-layer laminated 
cylindrical shells 
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Comparison of results for critical external pressure from the convex 
and numerical nonlinear programming 



Figure 6 Comparison of results for fundamental natural frequency from con- 
vex modelling and numerical nonlinear programming 
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Figure 7 Comprison of results for axial buckling load from convex modelling 
and numerical nonlinear programming 


curves are due to the change of vibration or buckling mode. The effect of un- 
certainty in elastic moduli on both the critical external pressure and the natural 
frequency varies greatly with the laminate configuration and the number of lay- 
ers which make up the laminated shell. For instance, for the 10-ply shells under 
the external pressure, the variability of the critical external pressure /3 goes from 
the mi n i mum 5.5% at 9 — 0° to the maximum 10% at 9 = 90°, whereas, for the 
5-ply shells under the same loading conditions, (3 varies between 5.9% at 9 = 0° 
and 11% at 9 = 47%. In comparison, the variability of the fundamental natural 
frequency due to the uncertainty in elastic moduli is smaller. For example, the 
uncertainty in elastic moduli brings about a maximum 5% and 6% variability of 
fundamental natural frequency aroung its nominal value for the 10-layer and 5-layer 
shells, respectively. However, the lamination angle corresponding to the maximum 
variability is different: 9 = 90° for the 10-layer shell and 9 = 43° for the 5-layer 
shell. In order to evaluate the accuracy of the above simple, analytic convex mod- 
eling analysis, we perform a numerical analysis of the critical external pressure and 
the natural frequency for a set of 10-ply laminated cylindrical shells by nonlinear 
progra mm in g . The results are displayed in Figures 5 and 6, which demonstrate a 
very good overall agreement between the results from the two different methods. 
Thus, once again, analytic convex modeling is proved very effective. It can also 
be seen from these two figures that when the lamination angle 9 is less than 45°, 
the results from convex modeling almost overlap with those from nonlinear pro- 
gramming. When 9 is greater than 45°, convex modeling predicts a slightly bigger 
variability than nonlinear programming. As far as the natural frequency and the 
critical external pressure are concerned, a design based on convex modeling results 
seems more conservative. However, it appears difficult to generalize this conclusion 
since it has been found that sometimes convex modeling may under-estimate the 
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variability of certain structural behavior. For instance, for the same 10-layer com- 
posite cylindrical shells as above, if the loading condition is axial compression, the 
buckling load results generated by nonlinear progra mmi ng can lie slightly outside 
the variation pattern predicted by convex modeling (Figure 7). 

It should be brought to the attention of the designer that, whenever one is 
confronted with imprecisely measured elastic moduli or with a collection of scat- 
tered data on elastic moduli, it is advisable to evaluate the potential variability 
of the structural properties, since such a variability may lead to the degradation 
in the load-carrying capability and dynamic behavior of the structure. Besides, it 
is worthy to note that, in many circumstances, while the traditional deterministic 
analysis based on using nominal values may lead to under-designing of the struc- 
ture, and the probabilistic analysis may not be performed due to the insufficiency 
of the available data (especially when an extremely low probability of failure is re- 
quired for the designed structure), the proposed non-stochastic, convex modeling, 
method may serve as a viable alternative to both deterministic and probabilistic 
methods. 

7 Conclusion 

Both analytical convex modeling and numerical nonlinear program min g are 
utilized to evaluate the variability of buckling load and natural frequency of the 
composite shells due to the uncertainties in elastic moduli. Based on the experi- 
mental data for the elastic moduli, the lower and upper bounds of the buckling load 
and natural frequency are evaluated. The bounds are very useful in practice and 
could be directly incorporated into design. It is remarkable that this simple non- 
stochastic, convex modeling enables the engineer to predict accurately the variation 
profile of certain structural behavior due to uncertainty in moduli, even when the 
experimental data are limited and hence the conventional probabilistic approach 
can not be utilized. 
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Abstract — This paper deals with determination of the best ellipsoidal model fitting the available 
limited experimental data. The problem is defined as that of finding the minimum volume ellipsoid 
containing all experimental data. A general transformation matrix for the rotation of AT-dimensional 
coordinate system is first obtained by the Gramm-Schmidt orthogonalization procedure. The use of 
this matrix makes it possible to search in all possible directions to find an ellipsoid with a minimum 
volume. The general procedure is illustrated by examples in which the real data is utilized. An 
invariance property of the response with uncertain parameters of different physical nature is also 
discussed. 
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Keywords — Convex modeling, Experimental data, Uncertainty analysis. 

1. INTRODUCTION 

The availability of uncertain, limited, information for the parameters either in a structure or 
in a excitation to which the structure is subjected, or in both, is often encountered in various 
branches of engineering; this is partially due to high cost of the measurements. For this case, Ben- 
Haim and Elishakoff [1,2] developed a novel approach, dubbed as convex modelling, to analyze 
vibration and buckling of beams, plates, and shells due to uncertain excitation or uncertain 
geometrical parameters. When the excitation is of the stochastic nature with some imbedded 
uncertain but nonstochastic parameters, the method of random vibration must be combined 
with convex modelling. These considerations led Elishakoff and Colombi [3] to propose a new, 
hybrid probabilistic and convex-theoretic approach to analyze dynamic response of structures. 
In the special case when the set of uncertain parameters is an ellipsoid, closed-form solutions 
were derived for the upper and lower bounds of the mean-square displacements of the structure. 
The direct comparison of probabilistic and convex analyses was performed by Elishakoff, Cai and 
Starnes [4]. 

Uncertainty modeling by methods alternative to probabilistic modeling was dealt in various 
contexts by Leitmann [5], Chernousko [6] and Schweppe [7]. Applications to structures include 
papers by Drenick [8], Shinozuka [9], Deodatis and Shinozuka [10], Lindberg [11], Koyluoglu, 
Cakmak and Nielsen [12]. The analysis of structures based on convex model for uncertain para- 
meters consists of two parts: one is the formulation of a deterministic objective function, namely, 
the stress, strain or displacement of the structure; the other is modeling of uncertain parameters, 

The study reported herein was supported by the NASA Langley Center. This support is gratefully appreciated. 
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which are represented as belonging to a convex set. In this paper, we focus on the second part, 
namely, determination of convex model. It is obvious that the response or the buckling load 
of the structure depend on the associated convex model. Hence, the problem of determining 
the best convex model for a limited available information about the uncertainty parameters be- 
comes of paramount importance. In the previous studies [1,2], for numerical convenience, the 
axes of the ellipsoids appearing in the convex models for the uncertain parameters were taken 
as directed along the coordinate axes. Thus, the results obtained based on such ellipsoidal set 
may be somewhat conservative for engineering design, since the volume of this ellipsoid may not 
possess the minimal property amongst all possible ellipsoids which can be constructed. In this 
study, we abandon this restrictive assumption. A general transformation matrix for rotation 
of TV-dimensional coordinate system is first obtained by the Gramm-Schmidt orthogonalization 
procedure. The use of this matrix makes it possible to search in all directions to find an ellipsoid 
with a minimum volume. The general procedure is illustrated by an example in which real data 
is utilized. 


2. 7V-DIMEN SION AL ELLIPSOIDAL CONVEX MODEL 

Assume that there are TV uncertain parameters a, (i = 1,2,..., TV) describing either in the 
structural properties or in the excitation. These parameters constitute an TV-dimensional para- 
meter space, namely, a T = . . . , a jv}. Assume that we have limited information on these 

parameters, represented by M experimental points, (r = 1,2, . . . , M) in this TV-dimensional 
space. The convex model assumes that ail these experimental points belong to an ellipsoid 

(a-a°) T G(a — a 0 ) < 1, (1) 

where a 0 is the state vector of the central point of the ellipsoid, and G is its characteristic matrix 


9 11 

<0 

t-* 

to 

9lN 

921 

922 

92N 

.9N1 

9N 2 • ■ 

■ 9nn . 


(2) 


which determines the size and the orientation of the ellipsoid. The matrix G is diagonal only 
when the axes of the ellipsoid are directed along the axes of the coordinates. Once G and ao are 
found, the ellipsoid is determined. The best ellipsoidal convex model in the class of the ellipsoids 
is identified with the one which contains all given experimental points but has a minimum volume. 
The possible steps to search for such an ellipsoid include the rotation of the coordinate system 
and construction of an ellipsoid whose axes are along the principle axes in the rotated coordinate 
system. If we can rotate the coordinate system in all possible directions, there must exist a 
direction along which the ellipsoid has a minimum volume. The proposed algorithm to achieve 
this goal is described as follows. 


2.1. Transformation Matrix for Rotation of Coordinate System » 

We first construct the transformation matrix for rotation of an TV-dimensional coordinate 
system. The new coordinates b are related with original ones a as follows 

b = Tjva, (3) 


where Tjv is a transformation matrix 


Tat = Tjv {fix, 62 , ■ . • , On- 1 ) , 


(4) 
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which represents a N x N square matrix and is dependent on N — 1 parameters, 01,0%, . . . , 0n-i, 
for the general iV-dimensional case; for example, for the two-dimensional case there is one pa- 
rameter 0 1 in T2, whereas for the three-dimensional case there are two parameters Q\ and 62 
in T 3 . 

The transformation matrix Tm can be constructed by any set of orthogonal vectors. One of 
the approaches to generate these orthogonal vectors is the Gramm-Schmidt orthogonalization 
procedure, which is briefly presented as follows: assume Vfc (k = 1, . . . , N) to be a set of linear 
independent vectors. We first normalize one of the vectors, say, Vi 


Uj 



( 5 ) 


Let a new vector W2 be defined as 


W 2 = V 2 - ciUx. 


( 6 ) 


We require orthogonality of W 2 to Ui 


u7w 2 = 0, 


(?) 


from which the constant ci is obtained 


ci =u 7 v 2 . 


( 8 ) 


Again, we normalize W 2 to yield 


U 2 = 


W 2 


V^wjw, 


( 9 ) 


Analogously, the general form of the fc th orthogonal vector and its normalized form can be 
obtained as follows: 


k - 1 

Wfc = Vfc — (\jJ Vfc) U, 7^ 0, k<N, 

2=1 
W fc 


( 10 ) 


Ufc = 


\/Wfc T Wfc 


The initial vectors Vfc can be chosen from any set of linear independent vectors. Here a set of 
Vfc are chosen as follows: 


' COS 0 \ 


' 0 ' 

sin 0\ cos 62 



sin 0\ sin 0 2 cos 63 


8i,k — 1 

. 

► ; Vfc = < 


sin 6 \ sin0 2 • - -sin0^/_ 2 cos 8 n-i 



. sin0!sin0 2 • ••sin0jv_2sin0jv_i , 


, 0 . 


O<0fc-i<|, (i = l, 2 ,...,AT; k = 2,3, ... ,N), 


where 6 i >m is the Kronecker delta, and the components in Vi are chosen to be the spherical 
coordinates in iV-dimensional space [14]. Obviously, the vector Vfc satisfies the following relation 


VjVfc = 1; 


(* = 1,2,..., AT). 


( 12 ) 
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By using equation ( 10 ), a set of orthogonal vectors Ufc (k = 1 , 2 , . . . , N) is derived as follows: 

Ui=V i; U fe = |°^- 2 |; k = 2 , 3 , . . . , N, ( 13 ) 

where Ok-2 is a vector with k —2 zero components and U& is a vector with N-k + 2 components 

-sin 0 fc_i 
cos 6k-i cos 9 k 


U k ={ 


> • 


(14) 


cos 6k - 1 sin 6k • • ■ sin#;v _ 2 cos 0 /y-i 
cos0fe_i sin^fc • • •sin0^r-2sin0jv-i ) 

Thus, the transformation matrix T;v for the rotation of the IV-dimensional coordinate system 
reads 

T n (61,62, ■ ■ ■ ,6 n - i ) = [Ui,U 2 ,.. . ,Ujv] , (15) 

which is an orthogonal matrix. Prom the general IV-dimensional transformation matrix, the 
specific cases can be obtained. For N = 2 

~ cos 61 — sin 6\ 

sin 6\ cos 6\ 


T 2 ( 0 i) 


(16) 


For N = 3 


^2(61,62) — 


0 

— sin 62 
cos 62 


cos 6\ — sin 6\ 

sin 61 cos 62 cos 6\ cos 62 
[ sin 61 sin 62 cos 6\ sin 6 2 
The coordinate system for T 2 (#i) is shown in Figure 1 . Note that the T 3(8 1,62) can be repre- 
sented as the product of two transformation matrices 

T 3 ( 0 1 , 0 2 )=T< 1) (62) Tf ( 0 j), 

where 

"10 0 " 

0 cos 62 - sin 62 

_0 sin 82 cos 82 
"cos 61 —sin 6\ 0 " 

sin 61 cos 6\ 0 

.0 0 1 . 

which represent two consecutive rotations around axis 1 and axis 3 ', respectively (see Figure 3 ). 
The rotational coordinate systems corresponding to T^fo) and Tg 2 ^(#i) are given in Figures 2 a 
and 2 b, respectively. 


T 3 1] (82) = 


T 3 2) (6i) 


(17) 


(18) 

(19) 

( 20 ) 



r 


Figure 1. Rotation of coordinate system associated with T 2 (#i). 
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2.2. Ellipsoid in Rotated Coordinate System 

By using transformation matrix Tjv given in equation (4), the given M points in the rotated 
coordinate system will have their new coordinates denoted by (r = 1,2, . . . ,M). To obtain 
the ellipsoid, let us first examine an IV-dimensional box of the form 


|b - b°| < d, 


(21) 


which contains all M points. The vector of semi-axes d T = {di, d 2 , . . . ,djv} and the vector of 
central points b 0T = {£>?, b%, ■ ■ ■ , b% } of the “box” in the rotated coordinate system are given by 


d k = \[mex{b^}-min{b^}], 
b °k = \ [max {&i r) } + min {^ r) }] , 


(r = 1, 2, . . . , M; fc = 1, 2, . . . , N). 


( 22 ) 


We now enclose this box by an ellipsoid 


*:= i 


9k 


(23) 


where gk are the semi-axes of the ellipsoid. There are infinite number of ellipsoids which contain 
the box given in equation (21). Clearly, the best choice is the one with the minimum volume. 
The volume of an IV-dimensional ellipsoid is given by 


N 


V e = C N JJs*:, 


(24) 


fc=i 


where Cn is a constant, depending on the dimensionality of the ellipsoid; for example, C 2 = it, 
C 3 — 4tt/ 3, etc; n denotes the product. The surface of the ellipsoid should pass through the 
corner points of the “box” (see equation (21). Therefore, 


N d 2 
fe= 1^ 


(25) 


We are interested in minimizing the volume V e of the ellipsoid, subject to constraint (25). We 
use the Lagrange multiplier technique. The Lagrangian reads 



N 

L = C N Y[g k 

fc= i 

+ a( 


(26) 

By requiring 

dL n 
dg k ~ ’ 

(* = 

1,2,... ,1V), 

(27) 

we obtain a set of equations 





Cn 

N d? 

n 9k ~ 

= 0, 

(i = 1,2,. -.,1V). 

(28) 


Multiplying equation (28) by gi and summing up the results with respect to i, we obtain 

N d 2 

NV,-2X^2%=0. 

i= 1 


(29) 
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Combining equations (25) and (29), we arrive at 


A = fv«. 

(30) 

Substitution of equation (30) into equation (28) results in 


— - NV e % =0. 
9i 9i 

(31) 

Since V e is nonzero, we get 


g i = y/ r Nd i , (i = 1,2,..., IV). 

(32) 


Thus, once the size of the box equation (21) is known, the semi-axes of the minimum-volume ellip- 
soid enclosing the box of the experimental data are readily determined by utilizing equation (32). 
If there are no experimental points at the corner of the box, the size of such an ellipsoid may 
further be reduced until one of the experimental points reaches the surface of the ellipsoid. The 
semi-axes of the ellipsoid in this case may be replaced by rjgk, where the factor rj is determined 
from the condition 


max 

r 

If there are some experimental points in the corner of the multidimensional box, the factor 77 
equals unity. The ellipsoid (23) can be rewritten in the form 

{b — b°} T D {b — b 0 } < 1, (34) 

in which b° is the vector of central points whose components are given by equation (22), and D 
is a diagonal matrix 

D = diag {(t? 0 i )~ 2 , {r}g 2 )~ 2 , . . . , (Wat)" 2 } • (35) 

2,3. The Ellipsoid with Minimum Volume 
The volume of the ellipsoid now reads 

N 

v e = C N T) N n 9k, (36) 

fc=i 

which is a function of a set of parameters 6k (& = 1,2, . . . , N — 1). Therefore, the best ellipsoid 
among these ellipsoids is the one which contains all given points and possesses the minimum 
volume, i.e., 

V e = min (V e (6 1 ,6 2 ,...,6 N - 1 )}. (37) 

01»02v»0JV — 1 

A possible approach to determine this ellipsoid is to search among all possible cases by increas- 
ing 6k (fc = 1, 2, . . . , N - 1) from 0 to tt/ 2 in sufficiently small increments A 6k, and to compare 
the volumes of so obtained ellipsoids. Once one finds the ellipsoid with minimum volume in one 
direction, say 6% (k = 1,2, N — 1), the ellipsoid can be transformed back into the original 
coordinate system by applying the transformation matrix Tjv- Hence, the vector a 0 of central 
point and the characteristic matrix G in equation (1) become 

a 0 = T^b°, G = T)vDTa/, (38) 

where Tjv = Tjv (0j, . . . , 6%_ 1 ) is given by equation (15); vector b° and matrix D are given by 

equations (22) and (35), respectively. It can be shown that G is a symmetric and positive-definite 
matrix which is nondiagonal in the general case. 


N 


HO 


£- 

fc= 1 


% 


.1 2 


9l 


<1; (r = 1,2, .. . ,M). 


(33) 
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3. NUMERICAL EXAMPLES AND INVARIANCE PROPERTY 

As an example, a set of uncertain parameters obtained from real tests on shell bucking [14] are 
chosen. These parameters represent Fourier coefficients of the half-wave cosine and half-wave sine 
representations, respectively, of the initial imperfections of shells. Chosen four parameters form 
the four-dimensional space. The values of Fourier coefficients of initial geometric imperfection 
derived in eight tests are represented by eight points in this space. They are listed in Tables 1 
and 3, respectively. 


Table 1. The values of uncertain parameters Af. for half-wave cosine representation. 


k 

i 

2 

3 

4 


1.800 x 10~ 2 

-5.000 x 10~ 3 

6.700 x 10~ 2 

7.000 x 10“ 3 

4 2) 

3.400 x 10“ 2 

-3.000 x IO -3 

0.653 

2.800 x 10- 2 

4 , 3) 

2.300 x 10- 2 

-6.000 x 10- 3 

8.300 x 10- 2 

2.000 x 10- 2 

4 4) 

1.100 x IO -2 

2.000 x 10- 3 

-2.300 x 10~ 2 

0.000 

4 5) 

2.000 x 10“ 3 

1.000 x IO -3 

1.600 x 10- 2 

0.000 

4 6) 

2.000 x 10- 3 

0.000 

2.400 x 10- 2 

0.000 

4 7) 

3.000 x IO -3 

0.000 

6.600 x 10~ 2 

1.000 x IO -3 


Table 2. Minimum volume of ellipsoid obtained by different increments, Ad. 


Increment Ad 

Orientation di 

Volume A e 

di 

02 

03 

5° 

90° 

90° 

o 

O 

iO 

1.7623 x 10“ 6 

3° 

90° 

3° 

87° 

5.6886 x 10~ 7 

2° 

90° 

2° 

88° 

4.4827 x 10- 7 

1° 

90° 

2° 

87° 

4.2200 x 10~ 7 


Table 3. The values of uncertain parameters for half-wave sine. 


k 

1 

2 

3 

4 

B l 3) 
b ! 4) 
b£ ] 

B? 

3.700 x 10- 2 

2.600 x 10- 2 

5.600 x 10~ 2 
2.900 x 10~ 2 

9.000 x 10“ 3 
-2.000 x IO" 3 

3.000 x 10~ 3 

-1.600 x 10" 2 
-1.000 x 10~ 2 
-1.000 x io- 2 

8.000 x 10 -3 

4.000 x 10~ 3 

5.000 x 10~ 3 
5.000 x 10" 3 

6.600 x IO" 2 
0.611 

7.500 x 10- 2 
-1.900 x IO" 2 

6.000 x IO" 3 

2.000 x 10" 2 
4.900 x IO" 2 

9.000 x 10- 3 
4.500 x IO - 2 

-6.000 x 10“ 3 

1.000 x IO" 3 

1.000 x IO" 3 
0.000 

8.000 x 10 -3 


Table 4. Minimum volume of ellipsoid obtained by different increments of Ad. 


Increment Ad 

Orientation di 

Volume A e 

0i 

02 

03 



5° 

25° 

0° 

85° 

5.2986 x 10" 6 

3° 

15° 

0° 

87° 

7.2538 x IO -6 

2° 

26° 

0° 

86° 

5.0534 x 10“ 6 

1° 

25° 

0° 

86° 

5.0479 x 10~ 6 
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Figure 4. A simply supported beam with concentrated loading. 


Tables 2 and 4, respectively, list the minimum volumes and orientations of the ellipsoid obtained 
by different increments A 6. The size of the increment is reduced so as to achieve the practical 
numerical convergence. The process of reducing A 6 is completed if the volumes of so constructed 
ellipsoids change by less than one percent. The minimum volume of an ellipsoid depends on the 
value of the increment A 6. The axes of the ellipsoid with the minimum volume are usually along 
the general, nonprincipal directions. 

An additional simple example of beam with simply supported ends and under two concentrated 
loads is shown in Figure 4. We are interested in the invariance property of convex modelling. 
To this end we assume that two hypothetical independent investigations use two different unit 
systems, for example, SI system and English system to analyze the same experimental data. The 
most and least favorable response of a structure based on the convex analysis read [2] 


^max 

Rmin 


| = R( a°) ± v / r T G -1 r, 


dR( a) dR( a) dR( a) 

dai ’ da 2 ’ ’ da^ 


a~a° 


(39) 


The invariance property of the convex model will assure that the responses obtained by two 
different unit systems are identical. Let us investigate a simple example. The bending moment 
at the midpoint of the simply supported beam (Figure 4) subjected to two concentrated loads 
Pi and P 2 reads 

M=±[P 1 L 1 + P 2 (L-L 2 )\. (40) 

First, we assume that two span lengths ai = L\ and a 2 = L 2 are uncertain parameters with a lim- 
ited measurements represented by four points in two-dimensional space as shown in Figures 5 and 
6 in different units, and P\ = 1 kN (0.2248 klb/), P 2 = 2kN (0.4496 klb/) and L = 3m (9.843ft). 
The ellipsoidal convex model and the maximum values of moment at the midpoint of beam are 
evaluated in two different units as follows: for SI system (m, kN) 


61° = 65°, a 0 = {1.0721m, 2.0003 m}, 
_ 37.1723 9.4686 ' 

G “ [ 9.4686 21.2821 J’ 


(41) 


and 

Amax = = 1.53575 + 0.27217 = 1.80792 kN • m, (42) 

for English system (ft, klb/) 


9° = 65°, a 0 = {3.5176 ft, 6.5631 ft}, 
[3.45585 0.88065' 

G — [0.88065 1.97794 ’ 


( 43 ) 
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Lx 

Figure 5. Ellipse of uncertain parameters L\ and L 2 in SI units (m). 


and 

77 max = = 1.13272 + 0.20072 = 1.33344 klb/ • ft. (44) 

Note that the unit factors for force and length A p and A l are 

Ap = 1 kN = 0.2248 klb/, 

(45) 

Ai = lm = 3.281 ft. 

Thus, 

X L XpM^H = 0.73757 x 1.80792 = 1.33346 klb/ • ft = (46) 

Is is shown that two results are identical, as expected. 

Consider now the case that two uncertain parameters have different dimensions. In order to 
maintain the invariance property of the responses with different units, nondimensional uncertain 
parameters are suggested to be used in the convex analysis. For example, assume in our pre- 
vious example that uncertain parameters are L\ = a\Xp and P 2 = o^Ap, where Oi and 02 are 
nondimensional uncertain parameters. The other, fixed parameters are L 2 = 2Ap, Pi = lAp and 
L = 3A l- Thus, the bending moment at the midpoint of the beam from (40) becomes, in view of 
equation (40) 


M (ai,a, 2 ) — - (ai 4- 02 ) ApAp = m (ai, 02 ) ApAl, 


F t (oi,a 2 ) = grad T {M ( 01 , 02 )} X P X L = f T (oi,o 2 )ApA L , 


(47) 





Figure 6. Ellipse of uncertain parameters L\ and L 2 in customary units (ft). 

where m(ai,a 2 ) and /(ai, 02 ) have nondimensional values. Since ai and a 2 are nondimensional, 
the a 0 and G (or G _1 ) obtained by convex analysis in equation (38) are also nondimensional. 
Hence, the maximum value of the response becomes, in view of equation (39) 

M max = M (a?, a°) + i/f t (a?, a§) 0-^(0?,^) 

f no / 1 (48) 

= lm (4,4) + y/ T (a?,a^)G- 1 /(a?,a§)| AiAp, 

where a 0 = {a?,*!®} ' s the nominal vector or central point of ellipse given by equation (38). As 
is seen the choice of the units can be arbitrary. 

4. CONCLUSION 

A general transformation matrix for rotation of JV-dimensional coordinate system was con- 
structed by using Gramm-Schmidt orthogonalization procedure. It was shown that the use of 
this matrix makes it possible to search in all directions to find the iV-dimensional ellipsoid with a 
minimum volume. Several numerical examples have been chosen for illustration. It is shown that 
the axes of the ellipsoid with the minimum volume have the general orientation in the parameter 
space. The invariance property of the response of structure with uncertain parameters of different 
units was investigated. It is shown that, with the nondimensional formulation, the invariance 
property of the response is retained. 
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Abstract— Buckling problems of two types of multi-span elastic plates with transverse stiffeners are 
considered using a method based on the finite difference calculus. The discreteness of the stiffeners is 
accounted for. It is found that the torsional rigidity of the stiffener plays an important role in the 
buckling mode pattern. When the torsional rigidity is properly adjusted, the stiffener can act as an 
isolator of deformation for the structure at buckling so that the deflection is only limited to a small 
area. Copyright © 1997 Elsevier Science Ltd 


1. INTRODUCTION 

Multi-span plates are used in many engineering applications where different parts of the 
plate are connected with one another and strengthened by stiffeners and interior supports 
to cover a large space. When the supports and stiffeners are equi-distantly spaced and 
all the constituent components look exactly alike so that the structure is periodic, the 
analytical finite difference calculus [1] is usually an appropriate method to use. However, 
the use of this method is normally l i m i ted to perfectly periodic structures. The finite 
difference calculus fails if the periodicity is disturbed due to misplacement in the location of 
the stiffener or support, which may find its way into the structure through the manufactur- 
ing process. 

In reality, although some periodic structures are designed to be completely identical for 
every constituent unit, they are seldom perfectly periodic. The deviation from complete 
periodicity is commonly known as disorder or irregularity. In this paper, we will address a 
particular type of disorder, namely, misplacement of stiffener or support. 

It has been found that in the presence of small misplacement of stiffeners or interior 
supports, the buckling mode shapes of certain structures display a localization in the 
neighborhood where an irregularity occurs. Localization phenomenon was first uncovered 
by the Nobel laureate P. W. Anderson [2] in physics. Its occurrence in structures has 
recently attracted much attention. Among others, Pierre and Plaut [3] considered the 
two-span column case. A more general case, the multi-span column, was recently discussed 
by Nayfeh and Hawwa [4] through the transfer matrix method and by the present authors 
[5] using a method based on the finite difference calculus. The deterministic buckling 
localization in cylindrical shells was investigated by El Naschie [6-12] and Hunt et al. 
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[13-16], whereas probabilistic aspects were studied by Xie [17]. Localization in a beam on 
elastic support was studied recently by Ariaratnam and Xie in ref. [23]. As for the plate, 
we used, in our previous study [18], a method based on a direct integration of the general 
governing differential equation. Such a method is effective, because it can account for the 
discrete nature of stiffeners so that the localization phenomena of the buckling mode can 
be discussed. However, a major drawback of the method is the large amount of com- 
putational effort involved, especially when the number of spans is large, because each span 
of the plate is considered separately and a large determinative matrix results. 

In this study, we combine the above method with the finite difference calculus to discuss 
(1) the general /V-span plate with transverse stiffeners and interior supports, which is 
structurally periodic except that one of the spans of the plate contains a disorder; (2) a 
piecewise periodic multi-span plate whose constituent segments are periodic themselves if 
the plate is partitioned into several segments. Unlike some other studies [19] which 
assumed a zero rigidity in torsion for the stiffeners, the contribution from the torsional 
rigidity is considered in the analysis. The entire treatment is theoretically exact and leads to 
a relatively simple formulation. It is shown in these two types of multi-span plates that 
even a single disorder could be responsible for the localized pattern of buckling modes and 
that the localization phenomenon can also occur in the piecewise periodic structures, 
though traditionally localization is discussed in the context of non-periodic structures. 

This study represents a generalization of two previous investigations: (a) on one hand 
ref. [18] which dealt with multispan columns; and (b) on the other ref. [21] which was 
devoted to passive control in columns. This study deals exclusively with elastic plates. 


2. PROBLEM FORMULATION 


The differential equation of the deflection surface of the plate subjected to a uniform 
compression P in the x-direction (Fig. 1) is 

3 4 w ' 


D r 


3 4 w + ^ 3 4 w 
. 3x 4 


+ +P^2L = 0 (1) 

3x 2 3 y 2 3/ / 3x 2 

where w is the transverse displacement, downward positive; D c is the flexural rigidity of 
the plate. For rectangular plates whose boundaries parallel to the x-axis are simply 
supported, the solution of eqn (1) can be represented in the following form [18] 


w’(x) = W(x)sin 


-(?)■ 


W(x) = A cos (ftx) + I? sin (ftx) + C cos (ftx) + D sin (ftx) 

(2) 

where A, B, C and D are unknown constants, which are to be determined by use of 
continuity and boundary conditions; and parameters ft and ft are denoted by 


ft = 
A = 


_7T 

b 
Pb 2 

7T 2 D r 


1 + 


IH) 


ft 


i 


i-)) 


( 3 ) 


here b stands for the width of the plate. 

Catering to a single, arbitrary ;th span, the solution can be written as 


W: 


”Hf)’ 


Wj = Aj cos (ftx,) + Bj sin (ft x,) + Cy cos (ftx,) + D ; sin (ftx,), 

0«Xft fly, (4) 
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Fig. 1. A multi-span continuous plate reinforced with transverse stiffeners. 


where is the length of the jth span and j ranges from one to N for an A-span plate. 

Here we consider a simplified case where there is a roller support under each interior 
stiffener. (The more complicated situation, namely, the plate without the interior support, 
is addressed in the Appendix.) In reality, we may deviate somehow from this condition. 
Instead of the presence of the interior supports, a more common occurrence in engineering 
practice is the use of girders or joists with plate structures, and sometimes the flexure 
rigidity of these girders can be so large that the deflection of the girders is negligible. If this 
happens, the deflection along the stiffeners can be regarded as zero. 

Using boundary conditions for an arbitrary y'th span 


( 1 ) 

( 2 ) 

(3) 

(4) 


Wj\x r o = 0 

w i\xr«i = ® 


dw, 


d Xj 
d Wj 


Oj-u Oj-i = ®/-i sin 


iry 


(5) 


d X: 




9j = © ; sin 


Try 


coefficients Aj, Bj, C y and D f can be determined with the aid of Mathematica [20] 


1 

A = —(lA sin (flidj) - A sin (/3 2 a y )]0 / + [-Acos(Aa ; )sin(Aa ; ) 

j i 

+ Acos(A«/)sin(A«y)]®/-i} 


Bj = ^-{P 2 [-cos (Picij) + cos (A«/)]®/ + [A cos (Aa ; ) cos (Aa ; ) - A 
+ A sin (A a /) sin (A a /)]©/-i} 
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Cj = —{[-ft sin (fttft + ft sin (ftft)]© ; + [ft cos (ft ft sin (/?!«/) 

$ 

- ft cos (ft ft) sin (ftft]© ; _i } (6) 


£). = i-{ft[cos(ftft - cos (ftft)]© ; + [-ft + ft cos(ftftcos(ftft 
S i 

+ ftsiftftftsiftftft]©^} 

ft = 2ftft[— 1 + cos (ft ft cos (ftft] + (ft; + /ft)sin(ft ft) sin (ft ft. 
The expression for the bending moment is 


M, 


~(B 2 W: d 2 Wj\ 

= -O e V — ± . 

\ dx 2 3y 2 / 


v ftr 

Using eqn (4), a moment-slope relationship can be established as follows 

D c, 


Mf_! = Mfto = m f-! sin -pS 

o 

M*" = M x \ x . =a . = m) sin^-; 


R 




L [Cl®;-l + C 2 ©y] 


m V = — [Cj©; + c 2 0;_i]. 


a, 


( 7 ) 


( 8 ) 


where Afjft, M*" are the bending moments at the two supports of the span, respectively. 
The superscript ‘R’ (‘L’) indicates that span of the plate is to the right (left) of the support 
in question (Fig. 2). The coefficients c x and c 2 are defined as 


/-I 



Fig. 2. Notations and positive directions. 
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Cl = —(fi\ - $)[- ft cos (ftft Sin (ft ft) + ft cos (ft ft sin (ft ft] 

S > (9) 

C 2 = - ^ 2 )[ft sin (ftft - ft sin (ft a ; )]. 

S i 


If a number of spans of the plate have the common length = a and are made of the 
same material, then the finite difference calculus may be applied in the discussion of that 
part of the plate. Equilibrium at a typical support, r, reads 


M? - M) = G7— , (10) 

3y 2 

where GJ is the torsional rigidity of the transverse stiffener. Substituting eqn (8) into the 
above, we have 


{2c l + k)G r + c 2 (0 r+1 + 0 r _i) = 0, (11) 

where k = 7T (2) aGJ/b 2 D c . 

Introducing the shifting operator E which is defined as E ft = 6 l+l , eqn (11) can be 
rewritten as 


[c 2 (E + E -1 ) + (2 cj + /c)]0, = 0. (12) 

Equation (12) is a second-order finite difference equation, whose solution is obtained by 


letting 


© r = e^. 

(13) 

Substitution into eqn (12) results in 


2c ! + k 

cosh (ft = . 

2c 2 

(14) 


Three different cases may arise and deserve separate considerations: 


Case 1 ~(2c i + k)/2c 2 ^ 1 
The solution of eqn (14) has the following form 


ft, 2 = ±cr; 

The attendant solution for 0, reads 


a = cosh 1 


2 c t + k 
2 c 2 


0; = A e ai + Be~ ai , 
where A and B are arbitrary constants. 


(15) 

(16) 


Case 2 — (2c x + k)/2c 2 =£ — 1 


The solution to eqn (14) takes the form 


0i,2 = ±(<* + ft); 


and the solution for 0, is 


a = cosh 1 


2c t + k 
2c 2 


j = V(-1) 


(17) 


0, = [74 cosh (cri) + B sinh (ai)\ cos (iri). 


( 18 ) 
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Case 3 -Is - (2d + k)/lc 2 =£ 1 
The solution to eqn (14) now becomes 


0i,2 = ± a; 


and, for & h we obtain 


a = cos 


2 d + & 

2c 2 . 


( 19 ) 


0, = A cos (ai) + Bsin(ai). (20) 

Introducing the following notations 

r ) = e ar , f 2 (a, r ) = cosh (ar) cos (nr), f 3 (a, r) = cos(arr) 

( 21 ) 

gi(ar, r) = e ar , g 2 (a, r ) = sinh (ar) cos (jrr), g 3 (a, r) = sin (arr) 

the three cases have a unified expression 

®r,i = AJla, r) + B^giia, r), i = 1, 2, 3, (22) 

where 0 r , corresponds to the three different cases when the subscript i varies from 1 to 3. 

We will consider two different kinds of continuous plates. The first kind of multi-span 
plate is structurally periodic except for a single disordered span that contains an 
imperfection. The second kind is a two-piecewise-periodic plate, which means that its first 
<?th spans and the rest of the N-q spans are periodic per se but they do not have the same 
periodicity. 

Suppose that we have an N-span plate (Fig. 3), which is periodic except that its 
( q + l)th span contains a span length imperfection which makes that span slightly longer or 
shorter than the other spans of the structure, i.e. 


Simple supports around the periphery 



4 t A II A""±”"A II t 'A A 


0 1 
y 


2 9-1 q 9+1 N-2 N-l N 


Fig. 3. An A-span continuous plate with a single disorder in the (q + l)th span. 
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a, = a for i = 1, . . q, q + 2, . . N; a q+x = a* (a A a*). (23) 

For a particular case a* = a, we recover the original, perfectly periodic structure. 

To facilitate the solution of the problem, we can treat the entire continuous plate as 
being composed of three segments. The first <j-span periodic plate constitutes segment I, 
the ( q + l)th span, namely, the disordered span, constitutes segment II and the last 
(N - q - 1) spans of periodic plate represent segment III. Assume first that both segments 
I and III themselves contain a large number of spans. For segments I and III, the finite 
difference calculus is applicable due to their structural periodicity. As to the disordered 
span, segment II, a separate consideration should be made. By following this procedure, 
we construct a solution composed of three parts with each part corresponding to a specific 
segment of the plate. Continuity conditions between those different segments are utilized in 
combination with boundary conditions at the ends of the plate to establish an eigenvalue 
problem. 

For the first q spans of periodic plate, we perform the finite difference calculus analysis 

e r = e\ = ©'sm-^-; (r = 0, 1, . . ., q), (24) 

b 

where the superscript denotes the sequence number of the segment in question; 0* takes 
on one of the three forms represented by eqns (16), (18) and (20), depending on the 
physical and geometrical conditions of the segment. 

For the disordered span, recalling eqn (8), we have 


D, 


M q = -r[ci6o + c 2 0?] 

a* 

< +1 = ~[eX + cXl 


(25) 


Or, in another form, 


0 o 


ii = a* CiMg + c 2 M q+x 

n -2 -2 

D c Ci - C 2 

a ii _ a* CiMg+i + c 2 M 
_ — — 

D c c\ - c\ 


(26) 


where c x and c 2 are obtained from the expressions for c x and c 2 by formally replacing in 
eqn (9) with a*. The treatment of the last N - q - 1 spans of periodic beam is similar to 
that of segment I, 

6 S = 0“ 9 _ i; (s = q + 1, q + 2, . . ., N). (27) 

Consider now a plate simply supported at its two ends (other boundary conditions can be 
treated in a similar manner). Then the boundary condition at the left end of the plate can 
be represented as 


Mi 


GJ X 


9 2 0p 

8y 2 


or ( Ci + /^©o + c 2 ©i = 0, 


(28) 


while the boundary condition at the right end of the plate reads 
-m l n = GJi 9 N ~r l 


dy* 


or (c x + k x )® l N- q - X + c 2 0^ 9 _ 2 = 0, 


(29) 
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where GJ X is the torsional rigidity of the transverse stiffeners at the two boundaries. In 
order to have a purely periodic structure, the stiffeners at the boundaries should have half 
the stiffness of those interior stiffeners, i.e. GJ X = GJ/2 or k x = k/2. This ‘half-stiffener’ 
concept has been used widely in ref. [1]. 

Conditions of continuity between the periodic spans and the disordered span of the 
beam, namely, between segment I and segment II, are 


-M\ + M* = GJ^L or (c, + A:)© 1 , + c 2 © 1 ,_ 1 + m f = 0 

9 y d c 

d\ = On or 


Q i a* /_ a r a L \ 

01 ‘ + = °- 


(30) 


Analogously, the continuity conditions between the second and the third segments are 


-M q+1 + M^ +1 = GJ — or (cj + A:)©™ + d©* 1 

3 y 2 


hi a l 
D 9+1 


0? = 0™ or ©o 11 + 


a* ( _ a r _ a l \ n 

+ = 0 . 


(31) 


Using the unified expression, equation (22), the rotation angles in the first and third 
segments can be expressed as 

©I = AJiu, r) + B lgl (a, r) (r = 0, 1, . . ., q ) 

TTT (32) 

©™ 9 -i = M fi(a, s - q - 1) 4- B 2 gi(a, s - q - 1) (s = q + 1, . . ., N). 

Substituting the above expressions in the boundary conditions (28) and (29) and the 
continuity conditions (30) and (31), we obtain six homogeneous algebraic equations, 

[(Cl + fcM*, 0) + C 2 f(a , 1)]A X + [(d + kiMa, 0) + c 2 gi (a, 1)]B X = 0 (33) 

[(ci + k)fi(a, q) + c 2 f{a, q - 1 )]A X + [(c x + k) gi {<x, q) + c lgi (a , q - l)\B l + m* = 0 

(34) 

f(a, q)A x + gi (a, q)B x - (- -W - Cz (— W +1 = 0 (35) 

c\ - c\\ a f Cl - c\\ a ) 

[(ci + k)fi(a, 0) + c 2 f(a, 1)]A 2 + [(cj + k) gi (a, 0) + c 2gi (a, 1)]S 2 - m\ +x = 0 (36) 


Ai + Bi + 


c 2 / « , C 1 ( a \ -L 

Ti id m 9 + "1 2 m « +1 

cf + c 2 \ a ) d - C 2 V a ) 


0 


(37) 


[(ci + k x )f(a, N - q - 1) + c 2 fi(a, N - q - 2)]A 2 + [(c x + k x ) gi (a, N - q - 1) 

+ c 2gi (a, N — q - 2)]B Z = 0, (38) 

where 


m. 


— r m q a 


D r 


m 


9+1 


m\ +1 a 

D, 


(39) 


Thus, we have six homogeneous algebraic equations, which can be expressed in a matrix 
form as follows 


[FWW<5} 6 xi = 0 


(40) 
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where [F(A)] is the coefficient matrix, and {<5} T = {A x , B x , m*, rn q+l , A 2 , B 2 }. 
Another kind of continuous plate is characterized by 


fa for i = 1, . . ., q 

[a* for i = q + 1, . . N 


(a =£ a*). 


(41) 


The schematic representation of such a structure is displayed in Fig. 4. For this problem, 
the first periodic segment which consists of the first q spans of the plate may fall into one 
of the three different cases, whereas the second periodic segment which is comprised of the 
remaining N-q spans may present itself in another three different cases. Therefore, there 
might be nine separate cases altogether, which makes the search for a solution rather 
complicated. 

Recalling eqn (21), a unified solution for this piecewise periodic plate can be written as 


= fA x fi(a u r) + B^ifau r ), i = 1, 2, 3; r = 0, . . ., q 

r \A 2 fj{a 2 , r) + B 2 gj{oc 2 , r), j = 1,2,3; r = q + 1, . . ., N' 


Using the boundary conditions (simple supports at two ends) 

(1) (Mq) 1 = (2) ~(M L N ) n = GjJ^± 

3y 2 3 y 2 

and the conditions of continuity 

(1) e\ = 0?; (2) (M ^) 1 = (Mq) 11 , 

the following homogeneous equations are established 


(42) 


(43) 


(44) 


Simple supports around the periphery 



t II t II t 


0 1 
y 


2 q - 1 


q + 1 N-2 N-l N 


Fig. 4. A piecewise periodic continuous plate. 
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[(*4 + c i)ft( a u 0) + c 2 fi(a u l)]Aj + [(&! + cOgXtt!, 0) + c 2 gi{a u l)]#! = 0 (45) 

[{k + c x )fj(cx 2 , N - q) + c 2 fj(oc 2 , N - q - 1 )]A 2 + [(£ + c^gjifr, N - q) 

+ c 2gj (a 2 , N -q- 1 )]B 2 = 0 (46) 
q)A x + gi (a i, q - 1)B X - fj{a 2 , 0 )A 2 - gj (a 2 , 0 )B 2 = 0 (47) 

[cif(a u q) + c 2 f(a u q - l)]Aj + [c x gl<x x , q - 1) + c 2g i(a u q - 1 )]B l 

-(^)[c^(a 2 , 0) + c 2 fj(a 2 , 1)]A 2 + ^j[c 1 g ; (ar 2 , 0) + c 2gj (a 2 , \)]B 2 = 0, (48) 

where the sub-indices i and j take on the value of 1, 2 or 3, depending upon which 
particular case the segments fall into. Again, the above four equations can be written in 
matrix form 


[F(X)] 4x4 {S} 4xl = 0, (49) 

where [F(A)] is the coefficient matrix, and {<5} T = {A u B x , A 2 , B 2 }. 

Both eqns (40) and (49) are homogeneous algebraic equations. Non-triviality of {<5} 
requires that the determinant of the coefficient matrix vanish 

Det[F(A)] = 0 (50) 

which constitutes a transcendental equation from which the non-dimensional buckling load 
parameter A can be solved in terms of other geometric and material properties of the 
structure in question. Once the buckling load parameter A is determined, we can use eqn 
(4) to calculate, span by span, the buckling mode shape for the entire structure, after the 
type of case has been ascertained. 


3. NUMERICAL EXAMPLES AND DISCUSSION 

In this section, we discuss the buckling load and mode shapes of the two different types 
of multi-span plates described in the last section. As the first example, consider the 
following case: 

a /7* 

— = 1, — = 1.1, N = 11, q = 5. 

b a 

As is shown in the above data, the plate consists of eleven spans, of which the sixth span 
contains a length imperfection which makes that span a bit longer than the other spans. 
Numerical results show that such an imperfection has a slight degrading effect on the 
buckling load. For instance, when k, the parameter characterizing the torsional rigidity of 
the stiffener, equals five the buckling load parameter A is 5.06. Compared with its 
counterpart of the periodic plate, which is A = 5.26, the reduction rate is only 4%. The 
buckling load reduction remains almost unchanged with the torsional rigidity of the 
stiffeners. Even when the torsional rigidity doubles, the reduction rate only amounts to 
5%. So, the buckling load decrease induced by the presence of the imperfection is not 
significant. However, the buckling modes are appreciably different for the plates with and 
without the imperfection (Figs 5-7). Moreover, as k increases, the buckling mode of the 
disordered plate becomes increasingly localized (Figs 6 and 7). The overall behavior of such 
plates is very similar to that of the continuous beams with torsional springs discussed in our 
previous paper [5], despite the difference of structural dimensionality between beams and 
plates. 
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0.06 -i 



Fig. 6. Buckling mode for a disordered 11-span plate ( k = 5). 


The second example is a 10-span, piecewise periodic plate whose first five spans have the 
length of a and the other five spans have the length of a* (assuming a* > a), and the plate 
is reinforced by transverse stiffeners with torsional rigidity of k = 10. When the ratio of a* 
to a equals unity, the plate reduces to a purely periodic plate. Here we discuss the case 
where the ratio is different from unity, say, a* /a = 1.1. As far as the buckling load is 
concerned, the difference between such a structure and the purely periodic plate is 
relatively minor. For the above structure, the buckling load parameter A equals 5.35, while, 
for the corresponding, exactly periodic plate, A is 5.77. So, the difference in the buckling 
load between the two is only 7%. More significant, however, is the difference in the 
buckling mode. For instance, when a*/a = 1.1 and k = 10, the deflection of the structure at 
buckling is largely confined to the left end, while those spans of the plate near the other 
end hardly experience any deformation (Fig. 8). Figure 9 shows the buckling mode of the 
same type of structure but with weaker stiffeners (k = 5). It is observed from this figure 
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Fig. 7. Buckling mode for a disordered 11-span plate (k = 10). 



Fig. 8. Buckling mode for a piecewise periodic plate (k — 5). 


that with a decreased torsional rigidity, the localization of buckling mode becomes less 
severe, and that deflection spreads from the right end to the left in a gradually attenuating 
fashion. 

From the examples shown above, it is well demonstrated that the torsional stiffness of 
the stiffeners should not be ignored in the investigation, just for the sake of simplification 
in analysis. As it turns out, the torsional stiffness plays quite an important role. It not only 
boosts the strength of the structure, but also localizes the loss of geometric rigidity of the 
structure at buckling to a small area so that any damage, should it occur, is kept to a 
minimum. 
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Fig. 9. Buckling mode for a piecewise periodic plate ( k = 10). 


For the plates contemplated here, small structural irregularities do not significantly alter 
the load-carrying capacity. However, the presence of such irregularities confines the 
buckling pattern associated with large deflection to a limited fraction of the structure. In 
this regard, the effect of such irregularities on the buckling loads can be considered 
favorable. As Nayfeh et al. [21] pointed out, by means of inducing ‘deliberately’ some 
irregularities in the system, one may confine the structural buckling to a limited part of the 
system only, which can be regarded as a passive control of the buckling process. 

The method used here seems to have a wide application, especially in continuous plates 
with a large number of spans. In fact, the larger the number of spans, the more 
advantageous the present method is over the other traditional methods such as the one 
based on the integration of the governing differential equation [18, 19], since it leads to a 
determinative matrix that can be orders of magnitude smaller than those necessary in other 
methods. Besides, the solutions generated by the present method are analytic and exact, 
and thus can be used as benchmarks for other numerical methods. 


4. CONCLUSION 

The method of investigation presented here can be applied to discuss essentially periodic 
structures where irregularities only occur in some local areas. It takes into consideration 
the discreteness of the stiffeners and, in particular, their torsional rigidity. As it turns out, 
the torsional rigidity of the stiffeners is important and should not be ignored in the 
discussion of the buckling mode shape. Adjusting the torsional rigidity of the stiffeners, one 
can achieve the goal of localizing the deflection of the structure at buckling to a small area. 
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APPENDIX 

This appendix describes how to analyze the buckling of periodically stiffened plate without the interior supports 
which were used in the main body of this paper. For such a plate, the deflection of the plate at various stations 
has to be accounted for. Accordingly, the coefficients A,, Bj, Cj and £>,■ in eqn (4) are now expressed, again using 
Mathematica [21], in terms of not only the rotation angles, but also the deflection 

Aj = -^{[- A A cos (An,) cos (An,) + Aftcos 2 (An,) - £2 sin (An,) sin (An,) 

+ Aftsin 2 (Aa,)]W/-i + [A A cos (An,) - AAcos(A aj)]Wj (Al) 

+ [A cos (A«,) Sin (An,) - A cos (/3, «,) sin (fta, )]©,_, 

+ [~A sin (A aj) + A sin (An,)]©;} 

A = 4-{[~ A A cos (A fl y) sin (A fl ,) + & cos (An,) sin (An,)] W h i 

+ [A A sin (A 0 /) - Pi sin (An,)] Wj + [A cos (An,) - Acos(A«;)]©/ (A2) 

+ [-A COS (An,) cos (A n,) + A cos 2 (An,) - A sin (An,) sin (An,) + Asin 2 (Aa,)]©;-i} 
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Cj = — {[0i0 2 cos 2 (0iay) - 0i0 2 cos (0ifly) cos ((ha,) + 0i0 2 sin 2 (0iay) - 0? sin (ft a,-) sin (0 2 fly)] W,_i 

+ [-ft ft cos (fray) + ft ft cos (ft a,)] Wj 
+ [— ft cos (0 2 fly) sin (ft a ; ) + ft cos (fta ; ) sin (0 2 fly)]0y-i 
+ [ftsin(ftay) - ft sin (ftfly)]0y} , 


(A3) 


where 


Dj = — {[0icos(0 2 fly)sin(0ia,) - ft ft cos (ft a,) sin (fta ; )] W}_i + [-)3fsin(ft« ; ) + ftft sin(fta ; )]l^ 
**/ 

+ [0i cos 2 (0,<jy) - ft cos (ft fly) cos (ftuy) + ft sin 2 (ft a ; ) - ft sin (fta y ) sin (0 2 a ; )]0y_, 

+ [-ftcos(ftfly) + ft COS (ft fly)]©;} 

Sy = 2ft ft - (0? + ^i) sin (ftay) sin (ftfly) - 2ftftcos(fta)cos (ftfly). 

The moments and shear forces are related to displacement and rotation angles as follows 

M?_ 1 =mjL 1 sin-H.; mf_, = ^ 

b a. 


(A4) 

(A5) 


IWy = my-isin^-; 
V* i = vf_ lS in^; 

Vy = Vy" sin-^-; 

' ' f> 


ft©,-! + ft©, + — (ftw;-! + ftw;) 

a i 

m) = - A ft© y + ft©y_, - — (ftw;. - ftWy_j) 
fly fly 

v f-l = + ft 0 i+l - — (ft WJ--1 - ftw,) 

fly L fly 

[ft©y - ft@y-i - —(ft IV, - ftWJ-i) 

fly l fly 


(A6) 


where 


ft = -^-(^l “ ftb{ftsin [(ft - ft)fly] + ft sin [(ft - 0 2 )a,] - ft sin [(ft + ft)a ; ] + ft sin [(ft + ft)fly]} 


ft = -g(~Pi + (Si)[ft (ftfly) - ft sin (ftfly)] 


(A7) 


ft = y{ft[ftft ~ ftft cos (ftfly) cos (fta) - ft sin (ftfly) sin (ftfly)] 

+ 01 [01 02 - ftft COS (ftfly) COS (ftfly) - Sin (ft fly) Sin (ft fly)] 

+ -^-v[ 2 ftft - 2 ft ft cos (ftfly) cos (ftfly) - (01 + /?!) sin (ftfly) sin (ftfly)]} 

ft = yft 02 (/?i - 0 2 )[cos ( 0 !fly) - cos ( 0 2 fly)] 

ft = ^ftft(/3? - 02 ) {ft sin [( 0 i - ft)fly] + ftsin[(0! - ft)ay] + ftsin[(ft + ft)fly] - 0 2 sin[(0, + ft)fl]} 

ft = yft 02 ( 0 i - 02 )[ 0 i sin ( 0 i fly) - 02 sin (ftfly)] 

ft = — f-{ 2 0?02 + 2/3 iy Sl + 0ift-p-(8 - 4o) - ( 0 ? + 0?ft + 0,0| + 0j + 20?-^- + 40i0 2 -g- 

2 ^ V “ COS [(ft - ft)fly] + (0t - 0?0 2 - 0i0i + 0 4 2 

201 -Jr - 4 ftft|- + 202^- - 0?fiv + 20102-—V - 0^v|cos[(0, + ft)«y]J. 


+ 2 ££-ri£v- 2ftft2 


Equations of equilibrium now become 


A/f - Mjr = GJ^ 

3y 2 

a 4 w r 


V? - V, = El - 


(A8) 


3y 4 


If the segment of the plate in consideration is periodic, then we have ay = a. On substituting the relations 
defined by eqn (A6), the equations of equilibrium have the form 
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— (W r -i - W r+1 ) - 2 d,@ r - 4(0,-! + 9 r+ i) = ^*0, 
a b 2 D c 

4(©, +1 - 0 r _O - M2d 5 W r - d 6 (W r+1 + W r _i)] = ^pL. Wr . 

a b 4 D c 

Using the shifting operator E, the above equations can be rewritten as 

— (E - E -1 )W r + [2 di + k + d 2 ( E -1 + E)]0, = 0 


ad 4 ( E - E -1 )0 r - [2d s + h - d 6 ( E -1 + E)]W r = 0, 


where 


k = 


GJTr 2 a 
b 2 D c ’ 


h = 


EIiAa 2 
b 4 D r 


Here El is the bending stiffness of the plate. 
Canceling ® r from eqn (A10) leads to 

d 4 ( E - E- 1 ) 


+ 


2 d 5 + h - d 6 { E" 1 + E) 


W r = 0. 


(A9) 


(A10) 


(All) 


(A12) 


L2d! + k + d 2 { E” 1 + E) d 4 (E - E -1 ) 

This is a fourth order difference equation. Assuming the solution is of the form W r — e^ r , we have the following 
governing equation 


d 4 sinh ($) 


2(2 d\ + k) + d 2 cosh(<|>) 


■ + 


4 d s + 2h - d 6 cosh (ip) 
d 4 sinh (<p) 


e<P' = 0 


(A13) 


A cosh 2 (<t>) + B cosh ((/>)+ C - 0, 
where A, B, and C are coefficients defined as 

A = d\ - d 6 d 2 , B = d 2 (4d 5 + 2h) - 4(44 + 2k), C = (44 + 2h){4d x + 2k) - d\ 
Solving eqn (A14), we have 


cosh (</>) 


■o 


-B ± V(S 2 - 4BC) 

a 12 = — 

a 2 ' 2 A 


The solution for w, takes different forms for the different situations: 

If B 2 -4AC»0 

Wi = wP + WP 

A 1 e r i' + A 2 e - t i‘, yi = cosh -1 (a:), if a\ s 1; 
wP = ■ [A 1 cosh(y 1 i) + A 2 sinh (yp)] cos (jri)> Yi = cosh -1 (ari), if cx^ s — 1; 
A!Cos(yii) + A 2 sin(yii), y x — cos -1 (o'!), if -1 =s 1; 


wP = 


(A14) 

(A15) 

(A16) 


(A17) 


A 3 e y 2 ‘ + A 3 e -I, 2 ! , y 2 = cosh -1 ( 0 : 2 ), if a 2 5* 1; 

[A 3 cosh(y 2 i) + A 3 sinh(y 2 i)]cos(7ri), y 2 = cosh -1 (a 2 ), if a 2 =s — 1; 

A 3 cos(y 2 i) + A 4 sin(y 2 /), y 2 = cos -1 (ar 2 ), if — 1 =£ a 2 =s 1. 

If B 2 - 4 AC < 0 

W t - A j cosh (pji) cos (p 2 i) + A 2 cosh (pi /) sin (p 2 i) + A 3 sinh (p, i) cos (p 2 i) + A 4 sinh (p 4 i) sin (p 2 i) , (A18) 

where 


Pi = sinh -1 


Pi = cos 


z\ + z\ - 1 + V ((*1 + Z 2 - l) 2 + 4zi) 


Zl + Z 2 + 1 - V((zi + zl- l) 2 + 4 zi) 


1/2 


1/2 


(A19) 


Since the difference equation (eqn (A9)) is of fourth order, it is necessary to have four boundary conditions, i.e. 
two at each end, fully to specify the problem. If the plate is simply supported at both ends, then the boundary 
conditions are: 


At the left end (i = 0): (1) w 0 = 0; (2) Mq = G7— 

3y 2 

(3) w N = 0; (4) M l n = -GJ f . 

3y 2 


(A20) 


At the right end (i = N): 


(A21) 
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Substituting the appropriate form of the solution w, into the four boundary conditions, we establish a set of four 
homogenous algebraic equations in terms of four unknowns A ly A 2 , A 2 and A 4 : 

Ai 

[/(A)] 4x4- ^ =0. (A22) 

A 4 

The non-triviality condition dictates that the determinant of the coefficient matrix [/(A)] should vanish 

det[/(A)] = 0 (A23) 

from which the buckling load parameter A can be evaluated as the lowest positive root of the above transcendental 
equation. The elements of the determinant are evaluated by using symbolic algebra [20]. Limiting cases of 
buckling are compared with the comprehensive monograph of ref. [22]. Once the buckling load is known, the 
buckling mode can be readily calculated following the procedure described in Section 2. In the case that the 
stiffeners are not uniformly distributed over the plate, partition of the structure may be carried out in order to 
take advantage of the finite difference calculus, which requires the use of continuity conditions between different 
segments. 




Pergamon 


Int. J. Solids Structures Vol. 34, No. 28, pp. 3755-3767, 1997 
© 1997 Elsevier Science Ltd 
All rights reserved. Printed in Great Britain 
0020-7683/97 $17.00 + .00 

PII : S0020-7683(96)00230-2 


s' ff-Jf _ 


EFFECT OF THE THICKNESS VARIATION AND 
INITIAL IMPERFECTION ON BUCKLING OF *7 7#* 

COMPOSITE CYLINDRICAL SHELLS : ASYMPTOTIC p fj 
ANALYSIS AND NUMERICAL RESULTS BY r < 
BOSOR4 AND PANDA2 


YI-WEI LI, ISAAC ELISHAKOFF 

Department of Mechanical Engineering, Florida Atlantic University, Boca Raton, 

FL 33431-0991, U.S.A. 

JAMES H. STARNES, JR 

NASA Langley Research Center, Hampton, VA 23664-5225, U.S.A. 

and 

DAVID BUSHNELL 

Lockheed Missiles and Space Company, Inc., Palo Alto, CA 94304, U.S.A. 

(Received 18 April 1994 ; in revised form 1 March 1996) 

Abstract — This study is an extension of a previous investigation of the combined effect of axi- 
symmetric thickness variation and axisymmetric initial geometric imperfection on buckling of 
isotropic shells under uniform axial compression. Here the anisotropic cylindrical shells are inves- 
tigated by means of Koiter’s energy criterion. An asymptotic formula is derived which can be 
used to determine the critical buckling load for composite shells with combined initial geometric 
imperfection and thickness variation. Results are compared with those obtained by the software 
packages BOSOR4 and PANDA2. © 1997 Elsevier Science Ltd. 


1. INTRODUCTION 

Due to various factors in the manufacturing process, thin cylindrical shells may exhibit 
variations in wall thickness. In spite of the fact that buckling of uniformly compressed 
cylindrical shells has been studied intensively for the past several decades, the influence of 
thickness variation on the buckling load has seldom been studied. In the previous research, 
we have investigated the effect of thickness variation on the axial buckling of otherwise 
perfect isotropic shells (Koiter et al., 1994a) and imperfect isotropic shells (Koiter et at., 
1994b). These studies resulted in a conclusion that, although the thickness variation pattern 
in the form of the classical axisymmetric buckling mode may have some deleterious effect 
on the load-bearing capacity, the most detrimental effect of thickness variation occurs when 
the wave number of the axisymmetric thickness variation pattern is twice that of the classical 
buckling mode. Asymptotic relationships between the buckling load reduction rate and the 
thickness variation parameter were established for isotropic shells of non-uniform thickness 
(Koiter et al., 1994a, 1994b). 

The present study aims at the combined effect of axisymmetric thickness variation and 
axisymmetric initial imperfection on the buckling behavior of composite shells. We approach 
this problem by using Koiter’s energy criterion of elastic stability (Koiter, 1945, 1966, 
1980). Here, we consider the small axisymmetric thickness variation, and as a first approxi- 
mation, only terms up to the first order of thickness variation parameter are retained. The 
final product of this discussion is again an asymptotic formula which relates the thickness 
variation parameter and initial imperfection amplitude to the buckling load of the structure. 
Therefore, this study is a direct generalization and extension of our former investigation 
(Koiter et al., 1994a, 1994b) to the anisotropic case. The asymptotic formula obtained 
herein encompasses the isotropic shell. Comparisons with results obtained by the computer 
codes BOSOR4 and PANDA2 are provided. 
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2. FORMULATION BY THE ENERGY CRITERION 


The nonlinear strain-displacement relations for cylindrical shells are 


du 1 ( dvrV 


£ x - 


8x 2\dx J ’ 
dv w 1 /<?w N2 


^ dy + R + 2\8y J’ 

dv du dw dw 
^ xy dx + dy + dx dy ’ 


d 2 w 

k x = 

dx 2 

d 2 w 

Ky= ~W 

0 3 2 w 
K ^ = ~ 2 ~fady 


( 1 ) 


where x and y are the axial and circumferential coordinates in the shell middle surface ; u 
and v are the shell displacements along axial and circumferential directions, and w is the 
radial displacement, positive outward ; s x , s y and y xy are strain components ; k x , k } , and K xy 
are middle surface curvatures of the shell ; R is the radius of the cylinder. 

Thickness variation of the laminated shell invariably exists due to imprecision involved 
in the fabrication process. Here we discuss the case that the thickness variation is axi- 
symmetric and of uniform configurational nature : each lamina has the same variational 
pattern : 


K{x ) = h 0Jc 


^1 — £COS 



h 0 , k H{x) (k= \ ~ K) 


( 2 ) 


where h k and h (yk are the thickness and the nominal thickness for the k - th layer, respectively ; 
£ and Pi are the non-dimensional parameters indicating the magnitude and wave number 
of the thickness variation, assumed to be the same for all the constituent layers ; K represents 
the total number of layers in the laminate. At first sight, the perfect homology of the 
thickness variation may appear as a restrictive assumption. If the constituent layers are 
produced by the same manufacturing process according to the same specification, one can 
not rule out the existence of similar deviations from uniform thickness. Such an assumption 
may shed some light to the question of thickness variation and lead to a tractable analysis. 
However, most shells are manufactured by being wound on a mandrel. The inner wall of 
the shell would probably be fiat and the outer wall would have all the thickness variation. 
In future, some numerical results will be reported for a variable thickness case where the 
inner surface is at a constant radius and all the thickness variation occurs on the outer 
surface. Here we assume that the middle surface of the shell with thickness variation only 
(no geometric imperfection) forms a perfect cylinder. 

With the model presented in eqn (2), elements of the stiffness matrices [A], [5] and [D] 
for the laminated shell with variable thickness become 


A,j — (Qij)k(hk — h k _ i) — H(x) ^ (Q<j)k(ho,k—ho,k-i) — H(x)a ri 

k = ! k=\ 

= ( QMhl-hl- X ) = \ t (QMhik-hi k _ ,) = [H{x)fb i} 

Z k = 1 Z *=1 


Dij = xl(Qu)k(U-hl_ i) = ^[H(x)f XiQMhlk-hl^) = [//(*)] V 


K 

£ 

k= 1 


0J= 1,2,6) (3) 

where b tj and d tj are elements of stiffness matrices for the corresponding uniform laminate 
with thickness h () ; QA s are the transformed reduced stiffnesses of the individual lamina and 
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n = U m +U b +Q. 

or, with use of the constitutive relations (8) and (9) : 

1 


( 13 ) 


n = 


A , , c 2 + 2A , 2 s x e y + 2A , 6 s x y xy + 2A 26 s y y xy + A 22 s 2 + A 66 y 2 xy + D x , k~ 


+ 2D l2 K x K y + 2D ]6 K X K„ + 2D 26 K y K xy + D 22 K y +D 66 K 2 -N 0 


dw dw 0 
dx dx 


dxdy. (14) 


Substitution of eqn (1) into the above formula leads to the energy expression in terms of 
displacements u,v,w: 


n 4 


^ii 


du 1 / <5w\ 2 

dx 2\dx J 


+ 2 A i: 


du 1 / dwA 2 
dx 2\dx J 


dv w 1 / dw\ 2 
dy R 2\dy 


+ 2 A u 


duUdw '' 2 


dx 2\dx 


du dv dw dw 
dy dx dx dy 


+2A 2 


-tAa 


dv w 1 / dw 
dy R 2y dy 


du dv dw dw 
<3y + dx + dx dy 


du dv dw dw 
dy + dx^ dx dy 


+ A 22 
d 2 w d 2 w 


dv w l/3vtA 2 
dy R + 2\dy 


(d 2 w\ 2 

+ D\i (-—-r) +2D l2 

\dx 2 J dx 2 dy 2 


A 8 2 w d 2 w d 2 w 8 2 w ( d 2 w x2 

+ 4-^16 "TT" TTY f4Z) 2 6 — T - "b D 22 ( ~ " 
dy 2 dxdy dy 2 oxdy \dy 


+ 4 D 6 


/d 2 w\ 2 
5 \dxdy ) 


, dw dwn 


dxdy. 


(15) 


In Koiter’s energy criterion of elastic stability, variations of energy are performed at the 
fundamental (pre-buckling) state. 

The second variation of the energy for buckling modes is 


1 n kr 
Pi\u] = 2 


. du\ 2 dufdv w 

^'l-' +2Al2 Tx\dy + R ] + 2A 


dx J 


du fdu 

16 3 .. 3 .. 


dv 

dx\dy ' dx 


(dv w\(du dv\ (dv 
+ 2Al6 + R){fy + Yx) + A22 [dy + 


+D 11 


d 2 w 

dx 2 


+ 2 D l2 


d 2 w d 2 w 


d 2 w d 2 w 

+ 4/^26 TTY'S - a b4D 66 
dy 2 dxdy 


dx 2 dy 2 
(d 2 w\ 2 


+ 4 D lt 


viA 2 / du dv 

~ + Rj +A66 \di + te 

v d 2 w 

dx 2 dxdy +I>2 ‘-\dy 


d 2 w\ 2 
2 


^dxdy J 


, dw 


dxdy. 


(16) 


We will discuss the effect of the most critical type of axisymmetric geometrical imper- 
fection w 0 (x) = — ph 0 cos (2 px/R) (Koiter, 1963; Tennyson et al., 1971) where h 0 is the 
nominal thickness of the shell, y is the non-dimensional parameter giving the amplitude of 
the imperfection, and p is the wave number of the axisymmetric classical buckling mode, 
which is given by (Tennyson et al., 1971) : 
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do not depend on the thickness. In the following, use will be made of the transformed 
stiffness matrices [A*], [5*] and [£>*], which are related to the matrices in (3) as follows : 

[A*] = [A ]~ \ [5*] = [B][A], [£>*] = [D]-[B\[A*m (4) 


Thus 


A * _ 

" H(x) 


a fy, Bfj = H(x)b*Dfj = [H(x )] 3 d% 


(5) 


where [a*], [/>*] and [d*] are counterparts, in the uniform laminate, of the transformed 
stiffness matrices [A*], [2?*] and [£>*] . They are given by 


[a*] = [a]" 1 , [/>*] = [b][a], [d*] = [d\-[b]{a*}[b]. 


( 6 ) 


We will deal with symmetric laminates, for which there is no coupling between bending 
and extension. Thus, we have 


Bn = o, M = 0 0,7 = 1,2,6). 


(7) 


The constitutive relations for the anisotropic laminate are (Vinson and Sierakowski, 
1986) 


[AV 


1 A lt 

A\2 

A 16 ' 

e x 


< 


A 12 

A 22 

A 26 

' By 

(8) 

.N xy , 


'^16 

A 26 

A(,6 1 

$xy- 


M x 


(Du 

D i2 

i>,6 ' 

K x 


My 

) = 

D \ 2 

£>22 

D 2 6 

• Ky 

(9) 

M xy \ 


\D, 6 

T>26 

dJ 

K X y 



where N„ N y and N xy are stress resultants, M x , M y and M xy are bending and twisting 
moments, acting on the mid-surface of a laminate. 

Membrane strain energy of a laminated cylindrical shell of length L is 


U m 


(N x s x + N y Sy + N xy y xy )AxAy. 


( 10 ) 


Bending strain energy reads 

1 

U h = - 


{M x k x + M y Ky + M xy K xy )dxdy. 


(ID 


For the shell under axial uniform end compression N 0 , potential energy of the applied 
load takes the form 


n = 


1 

2 


. dw 5w 0 \ 2 , , 

""lto + 77j 4v ^ 


( 12 ) 


where w 0 is the geometric initial imperfection. 
Thus, the total potential energy is 
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R 


4 V^i 


(17) 


We supplement the second variation with the additional bilinear term due to geometric 
initial imperfection 


PnWo,u] 


-N 0 


dw d w 0 
dx dx 


dxdy. 


(18) 


The third variation of the energy reads 


P 3 [u] = 


1 

2 

+ 2Ai 6 

+ 2A 26 


L \a ^ U ( ^ w ' 


+ A , 


du / dwY 

cbc^dy J 


dx\dx 

du dw dw 1 / dwA 2 / du dv 

dx dx dy + J ydy dx J J 


dv wA dw dw 1 / 3vtA 2 / du dv 

dy + Rjdxdy + 2\dy ) dx 


2 + ( dw ) 2 f^ U + W 


dx J V dy R 


(dv w\(dw\ 2 „ (du dv\dw dw) 

+ A22 {Ty + ^R jfej +2 ^ 66 (^ + ^j^^r 


(19) 


We now assume that the buckling modes of the shell with a uniform thickness remain a 
good approximation for the buckling modes of the shell with small thickness variations. 
We are at least ensured that the critical load obtained in this was is, by the energy principle, 
an upper bound for the actual critical buckling load. 

According to the study of Tennyson et al. (1971), the following expression for the 
buckling mode can be adopted for the laminated cylindrical shell with the aforementioned 
axisymmetric initial imperfection w 0 : 


, 2 px _ px ny 

w = ocos— — + C„cos— cos — 
K K K 


(20) 


where b and C„ are constants, n is the number of waves in the circumferential direction. If 
we recall the shell equilibrium equations in terms of displacements u, v and w (Vinson 
Sierakowski, 1986) 


j L \ i 

L\2 

L l3 \ 

L\2 

L 2 2 

! 

L 23 

\L\3 

L 23 

L 33 1 





' 0 ' 


u 


0 


V 

= - 

d 2 w 1 

1 

w. 


\ 

|«S 


( 21 ) 


where operators L tj are 

d 2 
L\ i 


d 2 „ d 2 d 2 

o-u ttt +2«i6 da 66 . 
dx 2 dxdy dy 2 


d 2 d 2 d 2 

Ln-a 16 ^+(a n +a 6 


d_ 

'Jy 


T13 — — ( fll2 ^ + «26 )> p 22 — a 66 T ~ + 2u 2 6 a _. a _. + a 22 


i! 

’ dx 2 


' dxdy 


dy 2 ’ 
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L 2 


7? 


a 26 dx + a 22 



a 22 , d 4 8 A 

— - -\- di i |- 4J 1 ( 

R 2 dx A 8x 3 8y 


+ 2(d l2 + 2d 66 ) 


d A „ , 8 4 

b 4d 2 6 

8x 2 8y 2 8x8y 3 


+ rf 22 — (22) 


we can obtain the expressions for u and v as follows 


a 12 , . 2px „ . px ny 

“-~2^; bsm ir +Q - c - sm -R a>s R 

px ny 

v = K n C n cos~” sin (23) 


where 


«(p 2 ai i g 22 + n 2 a 66 a 22 - a\ 2 p 2 - a 66 a x 2 p 2 ) 
(a, ^ +a 66 « 2 )(a 66 /n 22 « 2 ) -(a 12 +a 66 ) 2 p 2 n 2 


pa 6 6(p 2 ai 2 -n 2 a 22 ) 

(a, j/i 2 +a 66 n 2 )(a 66 p 2 +a 22 n 2 ) - (a l2 +a 66 ) 2 p 2 n 2 


(24) 


It should be mentioned that in deriving solution (23), we used an assumption in the 
studies of Hirano (1979) and Tasi (1966) that the coupling stiffnesses A 16 , A 26 , D l6 , and D 26 
are zero. They are identically zero for cross-ply laminates. When the laminate is composed 
of many layers, these coupling stiffnesses are small and can be neglected. 

In our previous numerical analysis of composite shells with axisymmetric thickness 
variations (Li et al., 1 994), we have shown that, in the absence of the geometric imperfection, 
the thickness variation with wave number being twice that of the classical buckling mode 
(p t = 2 p) has the most degrading effect of the buckling load. This result is also observed 
for isotropic shells (Koiter et al., 1994a, 1994b). Now we are interested in the combined effect 
of the most critical geometric imperfection and the most detrimental thickness variation on 
the load-carrying capacity. 

Substituting eqn (20) and (23) into the second and third variations, we obtain, after 
retaining only the first order terms in e : 


C 2 nL 

Pity] = ~ 7 ~r [d 22 n A + 2d 12 n 2 p 2 +4 d 66 n 2 p 2 + d ] ip A +a 22 (l +K n n) 2 R 2 

47? 3 


27? 3 


- N 0 p 2 7? 2 + 2a j 2 ( 1 + K n ri) Q„pR 2 + a l , Q 2 p 2 7? 2 + a 6 6 {K„p + n Q „ ) 2 7? 2 ] 

16<W‘ (l - 1 - jeW. Wl - le)-4AT,Ji V 


(25) 


P u [u 0 ,u\ 


4bh 0 N 0 p 2 pnL 
R 


(26) 


p sty] = ■ 


bC 2 nL 
87? 2 


a 12 


1 + y s V 2 +4^1 + -s )(1 + K n n)p 


+ a 22 \ l--s n 2 - 1 

2 ) a u 


1 + l + flnf l-xeV 2 +4fln/6„( 1 + qs )-4a 66 {\+\e)np{K n p + Q n n)\. (21 j 
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The energy expression to be considered is 

P 2 [m] + P x , [Wq, w] +P 3 [«]. 


(28) 


The equations for the initial post-buckling behavior are furnished by setting the partial 
derivatives of the energy expression (28) with respect to b and C„ equal to zero : 


b_ 

R 3 


1 6rf Uj p 4 l--e 


3 \ a\ 2 R 2 




4N 0 h 0 p 2 p C 2 


1 + \^jp 2 +4^1 + ^ej(l +K n n)p 2 


+ °”(' - + 1 + 


+ 4a u p 2 Q„[ 1 + \ e j “ 4 «6 6 (^1 + ^ ^j n P( K nP + <2„«) j = 0 


and 


“7 [d 22 n 4 +2d i2 n 2 p 2 +4d 66 n 2 p 2 +d x x p 4 + n 22 (l + K n ri) 2 R 2 
2R ' 


(29) 


- N 0 p 2 R 2 + 2a, 2 ( 1 + #„«) Q n pR 2 +a X{ Qlp 2 R 2 + a 6f> (K„p + nQ„) 2 R 2 ] 


bC, 

+ ^' 2 


~ 1 + T e V + 4 f 1 + + K » n )P 2 +«2zf 1 - xeV 


V 2 , 


+ 


2 2 
! 


1 + 5 £ ) + ai 2 ( 1 _ 5 e ) / ’ 2+4ail/?3 ^' , ( 1 + ^ e 


—4a 6 6 ( \ + ^^jnp(K n p+Q„ti)^ = 0. 


(30) 


With the solution C„ = 0 from eqn (30), eqn (29) yields 

4N 0 h 0 p 2 pR 2 


b = 


16d,,/( l-|g^- Q| a 2 ^ fl-^ej + a 22 i? 2 fl-^£ (-4 N 0 R 2 p 2 


(31) 


Bifurcation buckling with respect to the asymmetric mode with amplitude C„ occurs at 


b = - — [d 22 n 4 +2d X2 n 2 p 2 + 4d 66 n 2 p 2 + d u p 4 +a 22 (l +K„n) 2 R 2 

K 


- N 0 p 2 R 2 +2a l2 (\ + K n n)Q n pR 2 +a u Q 2 p 2 R 2 +a 66 {K„p + nQ„) 2 R 2 ] 

l a \2 


■ 1 + +4^1 + j(l +K„n)p 2 


+ a 22 j 1-^8 )« 2 


2 2 
a l2 n 


l + ^Va 12 fl-^eV+4d! 11 /> 3 e/l + ^£ 


-4a 66 ( 1 + ^jnp(K„p + Q n n) 


(32) 
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Equating the above two expressions for b, we obtain the equation for the critical buckling 
load N 0 

\6d u p^-^-^~(^--~ej+a 22 R 2 (l-^£ ) j-4N 0 R 2 p 2 

x [d 22 n 4 +2d 12 n 2 p 2 +4d 66 n 2 p 2 +d u p 4 +a 22 (l +K„n) 2 R 2 
- N 0 p 2 R 2 + 2a 12 (l+K„ri)Q n pR 2 +a l , Qlp 2 R 2 + a 66 (K„p + nQ„) 2 R 2 ] 

— 1 + 

+ 1 + ^ + a i2^1 - +4a, i P 3 Q„(l + ^ 

-4u 66 ^1 + ^e\np(K„p+Q„n) j = 0. (33) 


-2 N 0 h 0 p 2 pR 3 la 12 


hy 


2 +4j l + ^s)(l +K n ri)p 2 


+a 22 \ 1~ ~e )n 


In solving eqn (33), integer search must be performed with respect to the circumferential 
wave number n to arrive at the lowest value of N 0 . 

We define the non-dimensional critical load parameter X (sometimes referred to as 
knockdown factor in the literature) as 


X = 


N r , 


(34) 


where N d is the classical buckling load in the absence of both initial imperfection and 
thickness variation. N c , is given by (Vinson and Sierakowski, 1986) 

( L V C n C 22 C 3 3 +2C 12 C 2 3C,3 — Ci 3 C 22 C 2 3 Ci ! — C] 2 C 22 


N cl = min - 

\mn J 


where 


Ci i C 22 — C 12 


(35) 


mn\ 2 


L J 


R 


Cll — A 11 — +A 66 I — , C 22 — A 22 „ +^66 T 


R 


mn 


mn 


C 33 = C>n ( — ) +2 (£>i 2 + 2Z> 66 )( — j ^ j +D 22 y R j + 


mn\ n 


L \R 


A, 2 / mn 


C i2 ~ (A 12 +A 66 ){ j ) ( d 5 Cl3 ~ R ( £ ’ Cl3 ~ o\t>) + B ^\j} ' 


A 22 (n 


R \ R 


R 


3. RESULTS AND DISCUSSION 

Axial buckling loads can be determined from eqn (33) for composite cylindrical shells 
containing a small axisymmetric initial imperfection and a small axisymmetric thickness 
variation. For practical purposes, the results thus obtained should be considered conserva- 
tive, since the most detrimental case of geometric imperfection and thickness variation is 
investigated. However, since we ignored in our derivation the higher order terms in s, the 
results from the present study should not be deemed accurate for shells having large 
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thickness variation. Also, the results may not be conservative from a design point of view 
because we are limited here to axisymmetric variations. 

As a numerical example, we discuss shells made of carbon/epoxy laminae, whose elastic 
moduli are E x = 13.75 x 10 6 psi, E 2 = 1.03 x 10 6 psi, v 12 = 0.25, G u = 0.42 x 10 6 psi. The 
shell is 6 inches in radius and 30 inches in length and is composed of ten equally thick 
layers, each being 0.012 inch thick. The laminate configuration is [6/—0/d/ — d/Q\ sym , with 
the fiber angle 6 varying from 0° to 90°. 

Solving eqn (33) numerically for the critical load N 0 with integer search performed 
simultaneously with respect to the circumferential buckling wave number n, and then non- 
dimensionalizing the result according to (34), we obtain the critical buckling load factor 2 
for different cases of thickness variation parameter e and imperfection amplitude /i. The 
results are plotted in Figs 1 and 2. The results obtained here confirm numerically the 
previous first-order asymptotic formula 


2=1 — e (37) 

which holds only for the axisymmetric buckling cases for composite shells without initial 
imperfection. It is interesting to note that as long as the axisymmetric buckling mode 
dominates, the buckling load reduction factor 2 remains constant, irrespective of the 
change in the laminate construction. However, once the shell has an axisymmetric initial 
imperfection, the buckling mode becomes non-axisymmetric, and the buckling load 
reduction is strongly influenced by the stacking sequence of the laminae. Figure 3 depicts 
the results of the buckling load factor 2 for shells of different laminate profiles, such as 
[457—457457—45745°]^ and [167 — 16716 0 /— 16°/16°]^ m , together with the results for 
corresponding isotropic shells. It can be seen from this figure that the load-carrying capacity 
of composite shells is sensitive to thickness variation, and especially sensitive to initial 
geometric imperfection. The imperfection sensitivity is comparable to that of an isotropic 
shell. Although axisymmetric geometric initial imperfections cause most of the buckling 
load reduction, further degradation in the load-bearing capacity of the shell due to axi- 
symmetric thickness variations should not be overlooked. 



Fig. 1. Effect of thickness variation and imperfection on the buckling load (laminate configuration 
[6l — 6/6l — 6l8]„ m thickness variation parameter e = 0.1). 
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Fig. 2. Effect of thickness variation and imperfection on the buckling load ([8/ — 6/61 — 8/9], rm , 

e = 0.2). 





Fig. 3. Buckling load reductions in shells with different laminate configurations (Case 1 : isotropic; 
Case 2: [45/ -45/45/ -45/45]^; Case 3: [16/- 16/16/- 16/16L,J. 


In order to check the accuracy of eqn (33), we used BOSOR4 (Bushnell, 1974), a 
computer code for stress, buckling and vibration of shells of revolution, to generate a set 
of comparable data for the non-dimensional critical load parameter X. Since the classical 
buckling load has been used to non-dimensionalized the critical buckling load, it is necessary 
to check the results from eqn (35) with their counterparts from the numerical software so 
that a common basis can be established for the follow-up comparison of results for non- 
dimensional critical load X. For this purpose, software package PANDA2, with use of 
either the shallow shell or Sanders’ theories (Bushnell, 1987, 1996) was run in order to 
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compute the classical buckling load N d . Predictions are plotted in Fig. 4, together with 
those from eqn (35) and those from BOSOR4, which is based on Sanders’ equations. 

Figure 4 shows that the classical buckling loads from different sources agree quite well 
except in the range 53° < 9 < 80°. For this range eqn (35) and PANDA2 yield similar 
predictions with the shallow shell “switch” turned on in PANDA2. However, a significant 
discrepancy exists between predictions from shallow shell theory and Sander’s theory. For 
53° < 6 < 80° the shallow shell theory is significantly unconservative. A more refined theory 
is required for the accurate prediction of “classical” buckling load, N ch in this range of 6. 

In the BOSOR4 models of the axisymmetrically imperfect shells, half of the 30 inch 
length of the cylindrical shell is represented, with symmetry conditions imposed at x = 15 
inches. The 15 inch long BOSOR4 model is subdivided into six segments in order to get 
enough nodal points for a good convergence and in order to represent accurately the 
sinusoidal variation which is equal to the axial wavelength of the axisymmetric buckling 
mode of the perfect shell. The sixth segment, adjacent to the midlength plane of symmetry, 
is half as long. BOSOR4 can handle orthotropic walls with meridionally varying thickness. 
The shell wall in the BOSOR4 model has the same constitutive matrix as the 10-layer 
laminated shell with fiber angle 9 = 16°. (Note: in the BOSOR4 model the same off- 
diagonal “anisotropic” terms in the integrated constitutive law are assumed to be zero as 
is the case in the theory presented in this paper.) 

Figure 5 displays the results of the non-dimensional critical load parameter X obtained 
from the asymptotic formula (eqn 33) and from BOSOR4 for a 10-layer composite shell 
(laminate configuration: [ 1 6°/ — 1 6°/ 1 6°/ — 1 6°/ 1 6°]^ m ) which contains both the initial 
imperfection and the thickness variation. It can be seen from this figure that the asymptotic 
formula predicts the knockdown factor quite accurately. 

Finally, it is worth mentioning that as a special case if we let 


Eh 0 1 — v 

a \ 1 = a 22 = “ 7 > a \2 = Va l 1 > a 66 = X a l 1> #16 =a 26 = 0 

1 —v- ^ 


£7z 3 1 — v 

d\\ — d 2 2 = ^ 2(1 2 ^’ di2 = vdn ’ d(>6 = 3 d i6 = d 2 6 = 0 (38) 
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where E is the Young’s modulus, and v is the Poisson’s ratio, and furthermore, if we select 
the asymm etric mod e at the top of the Koiter’s semi-circle (Koiter, 1980), that is, let 
p = n = 3(1 — v 2 )R 2 /2h 0 ] 112 , eqn (33) reduces to its counterpart in the isotropic shell 

case, 





Buckling of composite cylindrical shells 


3767 


(1 -A)(l -e-X)- 3 ^ 3(l 2 y2) Xn(l + = 0. (39) 

Equation (39) is identical to eqn (21) in our previous work (Koiter et al., 1994b) if the 
small term s/6 is ignored compared with unity. 

Figure 6 shows the comparison of results in the isotropic shell case using Koiter’s semi- 
circle and those using integer search with respect to the circumferential wave number n. It 
is seen that the agreement is excellent. 

One should stress here that in order to obtain good correlation of test and theory, the 
buckling load of perfect shells with nonlinear bending prebuckling effects should be cal- 
culated when the effects of boundary conditions become significant, as in the case for short 
shells and shells of intermediate length, for which the “boundary layer” length, (rt) l/2 
comprises a significant fraction of the entire length of the shell. 

Acknowledgement — This research has been supported by the NASA Langley Research Center, through Grant 
NAG-1-1310. 
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Abstract 

The methodology of homology design is investigated for optimum design of advanced structures, for which the 
achievement of delicate tasks by the aid of active control system is demanded. The proposed formulation of 
homology design, based on the finite element sensitivity analysis, necessarily requires the specification of external 
loadings. The formulation to evaluate the worst case for homology design caused by uncertain fluctuation of 
loadings is presented by means of the convex model of uncertainty, in which uncertainty variables are assigned to 
discretized nodal forces and are confined within a conceivable convex hull given as a hyperellipse. The worst case of 
the distortion from objective homologous deformation is estimated by the Lagrange multiplier method searching the 
point to maximize the error index on the boundary of the convex hull. The validity of the proposed method is 
demonstrated in a numerical example using the eleven-bar truss structure. © 1998 Elsevier Science Ltd. All rights 
reserved. 


1. Introduction 

A great deal of attention is paid to structure/control 
integrated optimum design for advanced structures, 
such as space structures, adaptive structures and smart 
structures, which are required to fulfill complicated 
missions in high quality by the aid of active control 
systems (see Fig. 1) [1-3]. The precise control of geo- 
metrical properties of these structures is stringently 
demanded for high quality of performance and sophis- 
ticated structural design methodology considering the 
interaction with active control system and mechanism 
should be devised. Homology design can be a candi- 
date of such structural optimum design for advanced 
structures, utilizing the concept of homologous defor- 
mation , , in which a prescribed geometrical property at 
a part of the structure is maintained before, during 
and after the deformation. The quality of resolution of 
a huge and precise radio telescope was greatly 
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enhanced by realizing homologous deformation keep- 
ing the shape of the reflector structure parabolic in 
spite of the deformation caused by its own weight [4, 5]. 
The control cost of an active system to adjust the 
shape of a parabolic reflector is greatly reduced and 
the resolution quality is ensured easily. 

Two formulations based on the finite element 
method have been derived for homology design in sta- 
tic problems [6,7], However, the formulations were 
given in purely deterministic problems, where the 
external loadings must be specified. The advanced 
structures are made very flexible owing to lightweight 
requirement. Thus, the homologous deformation under 
specified loading appears to be easily disturbed by the 
uncertain fluctuation of loading. The sensitivity of 
homologous deformation with respect to such uncer- 
tainties should be estimated prior to the application of 
homology design. 

In this study, the formulation of homology design 
with an illustrative numerical example using the ele- 
ven-bar truss structure is given prior to the discussion 
on the uncertain loading. Use is made of the finite 
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element sensitivity analysis. In the numerical example, 
the uppermost deck subjected to uniformly distributed 
vertical load as nominal loading is kept both straight 
and horizontal, before and after the deformation by 
homology design. The methodology for the estimation 
of disturbed homologous deformation is developed 
based on the convex model of uncertain loadings [8, 9], 
in which uncertainty variables are assigned to discre- 
tized nodal forces and the existence domain of the 
uncertainty variables is confined within a convex hull 
of hyperellipse. The worst case of homology design is 
estimated on the boundary of the convex hull as the 
point that maximizes the error index. The latter is 
defined as the square of the Euclidean norm of displa- 
cements error from objective homologous deformation. 
The validity of the proposed method for the worst case 
estimation of the homology design is demonstrated in 
a numerical example using the eleven-bar truss struc- 
ture after homology design. 


2. Homology design under deterministic loading 

2.1. Formulation by finite element sensitivity analysis 

For the sake of specificity, let us consider the static 
deformation of the linear and elastic structure discre- 
tized by the finite elements. The formulation for hom- 
ology design under deterministic loading is 
summarized first and is governed by the following stiff- 
ness equation in matrix form after the incorporation of 
the geometrical boundary condition, 

[K]{u} = {f] (1) 

where [K\ is the stiffness matrix, {u} is the unknown 
displacement vector, {/} is the external nodal force vec- 
tor and the degrees of freedom of discretized structure 
are denoted by L. In this formulation, we need a trial 
design changed to satisfy the constraint of objective 
homologous deformation, which is given as follows: 

(HO)} = {0} (2) 


Eq. (2) consists of J equations with respect to nodal 
displacements and represents the constraint of the 
objective homologous deformation. The way of rep- 
resentation is not unique even for the same objective 
homologous deformation. The success of the objective 
homology design greatly depends on the manner of 
determining adequate trial design and represents the 
homologous constraint neatly. For the sake of simpli- 
city of discussion, we treat only linear homologous de- 
formation, the constraint for which is given by linear 
equations of nodal displacements, in this study. The 
homologous constraint equation in this case is rep- 
resented in matrix form, 

[C]{u} = {d} (3) 

where [C] is the constant constraint matrix of J x L 
and {d} is the constant vector of J components. 

If we find an adequate trial design a priori, the 
nodal displacements obtained by solving the stiffness 
equation, Eq. (1), satisfy the homologous constraint in 
Eq. (2) and there is no need to change the trial design. 
In general, however, it is hard to find the transparent 
homology design. Therefore we are interested in chan- 
ging the trial design such that the objective homolo- 
gous deformation is achieved under specified loading. 
For the design change, we judiciously select structural 
parameters p m and assign the design variables a m in 
the form of: 

Pm ^ P m {\ + a-m) (4) 

The upper bar means the value of current trial design 
hereafter. The effect of the design change caused by de- 
sign variables is linearly approximated and the change 
of nodal displacement vector is expressed in the follow- 
ing linear form: 

M 

{u} = [u\ + ^{i4)a m (5) 

m—\ 

where is the displacement sensitivity vector in the 
first-order obtained by the finite element sensitivity 
analysis [10]. Substituting the linearly approximated 
nodal displacement vector of Eq. (5) into the homolo- 
gous constraint in Eq. (3), we obtain the governing 
equation of design variables in the form: 

M 

£[C]{K> m = {r/}-[C ]{«,} (6) 

m= 1 

The governing equation of design variables can be 
rewritten in matrix form: 

[AM = {b} (7) 

where [A] is J x M rectangular matrix, {a} is the de- 
sign variable vector to be determined and {b} is a con- 
stant vector. In case that M is greater than J, the 
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governing equation, Eq. (7), has solutions which in 
general, are not unique. On the other hand, if J is 
greater than M, Eq. (7) scarcely has any solutions. The 
prediction about the solution existence based on the 
size of the coefficient matrix [A] mentioned above 
suggests that the increase of the number of design vari- 
ables enlarges the possibility of success in objective 
homology design. 

The solution to Eq. (7) is efficiently handled with the 
Moore-Penrose generalized inverse [11]. The general 
solution of Eq. (7) is given as follows when the con- 
dition of solution existence given by Eq. (9) is satisfied: 

{« } = [Ar{b} + {[I]-[A]-[A]){h} (8) 

([/]-MM-){h} = {0] (9) 

where [A]~ is the Moore-Penrose generalized inverse 
of [A], [7] identity matrix and {/?} arbitrary vector [11], 
The first and second terms of the right-hand side of 
Eq. (8) are called the particular and complementary sol- 
ution , respectively. Eq. (8) indicates that the way of 
changing current trial design to realize objective hom- 
ologous deformation in general is not unique, since 
arbitrariness creeps in the general solution of the de- 
sign variable vector through the arbitrary vector {h}. 
We employ only the particular solution to determine 
the design variable vector for design change. The sol- 
ution thus determined is the solution of least design 
change assuming that the Euclidean norm of design 
variable vector becomes a minimum owing to the 
property of the Moore-Penrose generalized inverse. 

The objective homologous deformation is not rea- 
lized strictly by the above determined design variables, 
since the governing equation, Eq. (7), is derived based 
on the linear approximation of nodal displacements 
change. The deficient approximation can be compen- 
sated by the iterative renewals of structure, in which 
the hanged structure according to the design variable 
solution is used as the trial design of the next renewal. 
The iterative renewals are stopped when the homolo- 
gous constraint equation is regarded to be satisfied, 
that is, the Euclidean norm of (H(«)} becomes suffi- 
ciently small. 




mm 


Fig. 3. Deformed structure of initial trial design. 

2.2. Numerical example using eleven-bar truss structure 

An illustrative example of homology design is 
demonstrated by using the eleven-bar truss structure il- 
lustrated in Fig. 2. The cross-sectional area of all the 
members is set equal to 100 mm 2 in the initial trial de- 
sign and the material property of the members is 
characterized by Young’s modulus, E = 70.0 GPa. The 
structure is subjected to a uniformly distributed verti- 
cal load of l.ON/mm as nominal loading on the 
uppermost deck. The deformed structure of the initial 
trial design is depicted in Fig. 3. The measure for dis- 
placement is put at the upper right of Fig. 3. The 
maximum vertical displacement turns out to be 
0.133 mm in the initial trial design. 

The objective homologous deformation is set so as 
to keep the uppermost deck straight and horizontal 
before and after the deformation. The homologous 
constraint equation is given in the form of 

vi -v m 

V 2 - v m 

V3 - V m 

v 4 v m 

where v, denotes vertical nodal displacement of node i 
and v m means the average of vertical nodal displace- 
ments on the uppermost deck, that is, 
v m = Ol +v 2 + v 3 + v 4 )/4. 

The realized homologous deformation by changing 
the cross-sectional area of all the members is illustrated 
in Fig. 4. The vertical displacement of the uppermost 
deck turns out to be 0.0516 mm. The manner for illus- 
tration of the deformation is the same as that of Fig. 3. 
In Fig. 4, the cross-sectional area of the member after 
homology design is indicated by the number attached 
to each member. The increase in total weight is 58.7% 
in this example. Three renewals of trial design are 


(10) 



Fig. 4. Realized homologous deformation and changed cross- 
sectional area. 
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needed to ensure the sufficiently small Euclidean norm 
of {H(«)j, i.e. less than 0.001 mm in this case. One of 
the transparent homology designs for this example is 
constructed so as to make the cross-sectional area of 
vertical and outer horizontal members extremely large. 
This method inevitably results in an extreme increase 
of total weight. On the other hand, we take the course 
of the least design change from the initial design which 
is chosen judiciously. 


3. Convex model of uncertain loading 

The legitimate question arises on how sensitive the 
homology design is to uncertain fluctuation of the 
loadings. The behavior of the structure is described by 
the following stiffness equation in the state of the hom- 
ology design: 

L K\{u} = {f} (11) 

where [X| corresponds to the stiffness matrix of the 
structure experiencing homology design at the specified 
nominal loading {/*}. The lower bar and asterisk 
denote the value after homology design and the value 
under nominal loading, respectively. The fluctuation of 
loading is represented by introducing N uncertainty 
variables s„ assigned to N nodal force components in 
the form of 

/„=/„*(!+£„) ( 12 ) 

The external nodal force vector {/} is expressed in the 
form of Eq. (13) by taking account of the uncertain 
fluctuation of nodal force components, 

{/} = {/*} + [F]{s} (13) 

where {s} is the uncertainty variable vector which con- 
sists of N uncertainty variables. The amplitude of vari- 
ation [E]{e} from the nominal loading is assumed to be 
unknown-but-bounded by a convex hull of hyperellipse 
with respect to uncertainty variables given in the form 
of 

<7 2 {e} t [1T]{e}<1 (14) 


[IV\ is a symmetric and positive-definite matrix that 
defines the shape of the convex hull. The magnification 
coefficient q governs the expanse of the convex hull. 
[W] and q are determined from the available infor- 
mation. 


4. Worst case estimation of homology design 

4.1. Formulation for worst case estimation 

In the context of convex analysis the problem is 
posed as follows: find the most harmful amplitude of 
variation to maximize the distortion from objective 
homologous deformation, which is called the worst 
case. The fluctuation of nodal displacement vector is 
expressed in the linear approximate form of Eq. (15) 
by the result of sensitivity analysis with respect to 
uncertainty variables for the structure governed by the 
stiffness equation, Eq. (11): 

N 

{«} = {«*} + £{#* (15) 

n=] 

{«*} is the displacement vector under nominal loading 
after homology design, i.e. objective homologous de- 
formation, and {ufi is the rate of change with respect 
to e„. Obviously the homologous constraint in Eq. (3) 
cannot be satisfied under the uncertain fluctuation of 
loading. The left-hand side of Eq. (3) after substitution 
of the fluctuating nodal displacement vector expressed 
by Eq. (15) is defined as the error vector from exact 
homologous deformation. The error vector {e} is writ- 
ten in the following form: 

N 

M = £[C]{^}e„ (16) 

n=l 

To evaluate the magnitude of the distortion from 
objective homologous deformation, we introduce a 
concept of error index, which can be defined by using 
the components of the error vector. In this study, we 
define the error index A 2 as the square of the 
Euclidean norm of the error vector. The index can be 
summarized as 


Table 1 

Estimated error index by convex analysis 


Error index, A 2 
([x 10 -3 mm 2 ) 

Uncertainty variables 





«i 

£ 2 

s 3 

£4 

2.82 

0.707 

0.000 

0.000 

-0.707 

2.82 

0.500 


-0.500 

0.500 

1.94 

0.000 

-0.707 

-0.707 

0.000 

0.00 

0.500 

0.500 

0.500 

0.500 
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A 2 = {e} T [Z>]{e} * (17) 

where [D] is an N x N symmetric matrix defined by 
the rate of change of nodal displacements with respect 
to uncertainty variables. 

The worst case of uncertainty variables, which gives 
the maximum value of the error index, is searched 
within the convex hull of Eq. (14). As the error index 
is expressed in second-order form with respect to 
uncertainty variables and becomes the convex function, 
the point of worst case is located on the boundary of 
convex hull. The search is carried out by the Lagrange 
multiplier method employing the following functional: 

n = A 2 -% 2 {e} T [lE]{e}- 1) (18) 

where X is the Lagrange multiplier. The uncertainty 
variables in the worst case are so determined as to 
satisfy the stationary conditions derived as follows: 

^ = 2([D] - = {0} (19) 

^ = 9 2 { £ } T [fT]{ £ }- 1=0 (20) 

The stationary condition of Eq. (19) results in the 
eigenvalue problem formed by the eigenvalue Xq 2 , 
which is rewritten by p, and the eigenvector { £ }. On 
the other hand, the stationary condition of Eq. (20) 
gives the normalizing condition for the eigenvector. 
Generally, N eigenpairs are obtained by solving the 
eigenvalue problem of Eq. (19) under the normalizing 
condition of Eq. (20). In view of Eqs. (19) and (20), 
the error index of Eq. (17) becomes 

A 2 = { £ } t [T>]{ £ } = V( £ ) t [IF ]{ £ } =X = p/q 2 (21) 

as the magnification coefficient q is constant and the 
maximum value of the error index, i.e. the worst case, 
is estimated by the maximum eigenvalue obtained by 
solving Eq. (19). The eigenvector corresponding to 
maximum eigenvalue gives the uncertainty variable 
vector for the worst case. 




4.2. Numerical example to estimate worst case of 
homology design 

The eleven-bar truss structure obtained from the 
homology design in Section 2.2 is employed in this nu- 
merical example to evaluate the worst case of the dis- 
tortion from objective homologous deformation caused 
by uncertain fluctuation of the uniformly distributed 
vertical load. Four uncertainty variables are assigned 
to vertical nodal force components for nodes 1, 2, 3 
and 4 in Fig. 2. The fluctuation of the uncertainty vari- 
ables is bounded by hypersphere with radius 1.0, 
which is expressed by setting q equal to 1.0 and [W] as 
an identity matrix in Eq. (14). The nodal force com- 
ponents fluctuate from 0 to twice the nominal values 
in this case. 

The result of the worst case estimation obtained by 
solving the eigenvalue problem of Eq. (19) is shown in 
Table 1. The error indices estimated by eigenvalue, i.e. 
A 2 = j a/q 2 , are listed in the leftmost column and corre- 
sponding uncertainty variables, i.e. eigenvector com- 
ponents, are listed in the right columns. The error 
index is maximized by two cases of uncertainty vari- 
ables namely ( £ |, « 2 > £ 3 , £ 4 ) = (0.707, 0.0, 0.0, —0.707), 
called worst case 1 and ( £] , e 2 , £ 3 , e 4 ) = (0.5, —0.5, 
-0.5, 0.5), called worst case 2. The actual error indices 
evaluated by the analyses employing the nodal force 
components corresponding to the worst cases 1 and 2 
turn out to be equal to 2.82 x 10~ 3 mm 2 in both cases. 
The deformations of the eleven-bar truss in the two 
worst cases are illustrated in Figs. 5 and 6, respect- 
ively. The manner of illustration is the same as that 
employed in Fig. 3. 


5. Concluding remarks 

A formulation of homology design based on the 
finite element analysis is reviewed. The change of dis- 
placements by design variables is approximated in the 
first-order by means of the finite element sensitivity 
analysis and the governing equation for design vari- 
ables is derived in the matrix form with a rectangular 
coefficient matrix. The equation is handled with the 
Moore-Penrose generalized inverse to obtain the least 
design change solution. The illustrative numerical 
example using the eleven-bar truss structure demon- 
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strates that the uppermost deck is maintained straight 
and horizontal even after the deformation caused by 
uniformly distributed vertical load. 

A methodology to estimate the worst case of hom- 
ologous deformation caused by uncertain fluctuation 
of loadings is proposed by means of the convex analy- 
sis. The uncertainty variables are assigned to discre- 
tized nodal force components and the fluctuation of 
the uncertainty variables is bounded within a convex 
hull. The worst case to maximize the error index of 
homologous deformation is searched on the boundary 
of the convex hull. The search is carried out by the 
Lagrange multiplier method. The numerical example 
using the eleven-bar truss structure proves the validity 
of the proposed method by demonstrating the identifi- 
cation of the worst case cased by uncertainty variables 
confined within a hypersphere. 
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BUCKLING OF STRUCTURES WITH UNCERTAIN BUT NOT-NECESSARILY 
RANDOM IMPERFECTIONS - PERSONAL PERSPECTIVE 
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Dept, of Mechanical Engineering 
Florida Atlantic University 
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ABSTRACT 

The previous review on stochastic buckling of structures was 
written by Amazigo in 1976. This review summarizes some of the 
developments which took place in recent two decades. A brief 
overview is given of the effect on uncertainty in the initial 
geometric imperfections, elastic moduli, applied forces, and 
thickness variation. For the benefit of the thinking reader, the 
review has a critical nature. 

It should be noted that this manuscript has yet to be 
completed. The list of topics to be included is listed below. The 
final version of the manuscript will be obtainable by informing 
this writer. 


1. INTRODUCTION 

The general theory of buckling and postbuckling behavior of 
elastic structures was worked out by Koiter (1945, 1963). Further 
contributions were provided by Budiansky and Hutchinson (1964) , 
Arbocz (1985) and other investigators. For a bibliography the 
reader may consult, for example, with the articles by Budiansky 
(1974), Budiansky and Hutchinson (1979), Koiter (1985), and Arbocz 
(1990). 

There are many other investigations dealing with the chism 
that exists between the theoretical analyses and the experimental 
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results. Most unfortunately, the experimental results “misbehave” 
their way and do not match the theoretical predictions. In these 
circumstances it was not unnatural to look for the uncertainty as 
a responsible factor for the scatter in experimental results. One 
conceptually understands that there are not two identical shells 
produced by the same manufacturing procedure. Motivated by this 
idea, the investigators could ascribe the scatter in buckling loads 
to the scatter in initial imperfections. 

Next step made was to identify uncertainty with randomness and 
to utilize the probabilistic methods. We find the first hints of 
these thoughts in the paper by Hoff (1949) : 

"...The nature and the magnitude of the disturbance must be 
established from a statistical investigation of the conditions 
under which the structural element or part of machinery will 
be used. The safe of the system can be safeguarded if it is 
made stable for all disturbances which have a probability 
greater than a required minimum." 

This idea, apparently independently, was pursued by Bolotin 
(1958) . He postulated, in brief, that the buckling load X of a 
structure can be expressed as a deterministic function of a finite 
number of parameters representing the initial imperfections: 

X = <p{Z\t £ 2 # • • • / £jv) 

where N is the number of terms taken in expansions, we also assume 
that we are given a particular function <p, and that join 
probability density 
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( 2 ) 


fx(*l • ? 2 ' ‘ * * • %n) = 


n 


n (§* * x i * ?,• + a?,) 

/-i 


of random initial imperfection vector, denoted by 
X = (X 1 ,X 2 , . . . ,Xff) T i.e. the probability that the random components 


X { of the vector X will belong to the interval + A£,-) , whereA£ f 

is an increment. Note that the random variables are denoted by 
capitals, whereas the possible values they can take on are 
identified by lower-case notation. Bolotin applied this method to 
a cylindrical panel under uniform compressive load along its curved 
edges, with the initial imperfections represented by a single 
normally distributed amplitude parameter. A single-term Galerkin 
approximation yielded an equation of type (1) . Conceptually such 
an one-term analysis is not complicated. Once relation of type (1) 
is obtained, and the probability density of the initial 
imperfection X ^ is specified or assumed (i 0 is the index of the 


governing initial imperfection parameter) one calculates the 
reliability of the structure. We first note that due to assumed 
randomness of the initial imperfection, the associated buckling 
load turns out to be also a random variable, denoted by A. The 
reliability at the load level a is defined as the probability that 
the structure will not buckle prior to a, or in other words, it 
will live beyond "age" a : 

R( a) = Prob ( A > a) C 3 ) 

Having determined the reliability of the structure, one proceeds 
with its design as follows. One should have a codified reliability 
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r, i.e. the level of reliability below which the performance of the 
structure is defined as an unacceptable one. The probabilistic 
design criterion demands that 

R( a) > r 

Combining Egs. (3) and (4) leads to 

R( a) - Prob(A > a) > r (S) 

Inequality (5) leads to a possibility to solve some basic problems 
of stochastic buckling. If the left and right sides of Eq. (4) are 
known then one can check if the probabilistic design criterion (4) 
is met, or it is violated. If some probabilistic characteristic of 

the initial imperfection, say its variance d^, is unspecified, one 

can calculate its maximum admissible level of it max d^ , such that 
the design criterion is satisfied. The value max d^ is obtained 
by solving an equation 

R( a) =Projb(A>a) =r ( g ) 

Solution of this type of problems may then be introduced in the 
quality control measure; if the variance of the initial 
imperfection exceeds max v f the structure is declared unacceptable. 

This third problem consists in determining the design load a r , such 

that if a < a r then the reliability will not be less than r. 

The reliability of the structure at the nondimens ional load level 
a can be rewritten as 
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(7) 


R ( «) = Prob (- < X fo * fj) 


where f 1 is the value at which the limit load equals a. This 
implies that satisfies an equation: 

X = *>(! i) = a (8) 

Hence 

?i = <P~ l (ot) (9) 

If, for example X = X/ Q is a random variable having a normal 
distribution with zero mean and mean-square deviation d 


f*(x) = 




Then the reliability becomes: 


( 10 ) 


i?(a) = Prob(-<p~ l (a) <■ X < (a)) 


( 11 ) 


= 2 erf 


^ x (a) 



( 12 ) 


This allows to find the probabilistic design load a r , such that if 
a = a r than the least reliability of the structure equals r: 


a r 



(13) 


In order to illustrate the stochastic imperfection sensitivity 
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concepts let us consider a simple structure, namely, a column on a 
nonlinear elastic foundation 


El 


dV 

dar 4 


+ P 


dV 

dx 2 


+ - k 3 w 


3 = -i 


dV 

dx 2 


(14) 


where w(x) = initial imperfection, w(x) - additional deflection P 
= axial load, and k 3 - nonlinear spring coefficients of the 

foundation. The buckling of the perfect column on a linear 
foundation is a classical, textbook problem. The imperfect column 
on a nonlinear softening elastic foundation exhibits imperfections 
sensitivity in that the limit load the structure may support may 
turn out to be far less than that of perfect linear counterpart. 
Application of the Galerkin method for the column that is simply 
supported at its ends yields in a single-term approximation the 
following equation, derived by Fraser (1965) in his Ph.D. 
dissertation: 


(1 - X)» - ".E* 


32 


(15) 


where % m = initial imperfection amplitude associated with m half- 


sine waves in axial direction, X = nondimens ional limit load, s is 
a value depending on the physical parameters of the system. 

This type of analysis can be demonstrated on the imperfection 
sensitivity of a shell with a non-axisymraetric periodic 
imperfections, studied by Koiter (1963) : 
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( 16 ) 


where x = axial coordinate, y circumferential coordinate, g ~ non- 
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dimensional initial imperfection amplitude, i c = the number of 

half-waves at which the associated perfect shell buckles, L = shell 
length, h — shell thickness. Koiter (1963) arrived at the 
following equation relating the buckling load with the initial 
imperfection amplitude: 

(1 - X) 2 + 6 cg\ = 0 < 17 > 
when X = P tim /P c , where P tim is the limit load, P c = classical 

buckling load of the perfect shell. 

The analysis which is based on a single term approximation is 
quite similar in its general spirit but not in its particulars to 
the asymptotic analysis developed by Koiter (1945, 1963) and 
Budiansky and Hutchinson (1964) . The asymptotic equations or the 
equations based on a single-term Galerkin approximations can be 
utilized for explanatory purposes. 

Indeed, we are looking for highly reliable performance, 
associated with the probability of failure, say 10" 6 or even less. 
Realizing this, one immediately should cast a doubt on the 
possibility that the highly simplified expressions (which are of 
extreme importance to capture the physical phenomenon itself) would 
reliably produce the required high reliability. In other words, 
simplified expressions may be unusable to calculate extremely small 
probability of failure. 

It was perfectly valid to utilize a single-term Galerkin 
approximation in Bolotin's (1958) early work. Analogously, 
application of Koiter' s (1963) asymptotic expressions by Thompson 
(1967) , Roorda (1972) and Hansen and Roorda (1973) served a purpose 
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of illustrating the reliability approach in imperfection-sensitive 
structures. Yet it appears that an industrial firm, for example, 
cannot use simplified expressions to justify then reafter 
reliability calculations with attendant extremely small 
probabilities of failure; taking into account additional terms in 
Galerkin expansion or additional terms in asymptotic expansions, 
may significantly alter the resulting probabilities of failure, and 
invalidate the proposed designs. Yet, some very recent works still 
utilize the deterministic asymptotic expansions for reliability 
calculations (see, e.g., Cederbaum et al 1996). 

We must assume that these fine points were perfectly 
understood by some investigators quite early, since they did not 
follow the single-term or deterministic asymptotic methodologies, 
superimposed with treating the imperfection amplitude as a random 
variable. 


2. STUDIES BASED ON ERGODICITY ASSUMPTION 


In his review paper Amazigo (1976) stresses, relating to the 
equation (1) postulated by Bolotin (1958) : 

"It is however a nontrivial problem to obtain on (1) and 
perform the above analysis for n > 2, say. It is this 
difficulty that limits the effectiveness of this method." 
Instead of utilizing the concept of the random variable, as in 
works by Bolotin (1958) and Thompson (1967) , the scholars of the 
Harvard group correctly decided to adopt the theory for random 
functions, identifying the initial imperfections as random fields 
with specified probabilistic characteristics, namely, the mean 
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initial imperfection function and the covariance function. 
Apparently first such study 'dealing with imperfection sensitive 
structures was undertaken by Frazer and Budiansky (1969) . They 
studied the imperfect column on a nonlinear elastic foundation. 
The length of the column was taken to be infinity. The following 
assumptions were made about the initial imperfection field : (a) 
they were considered to form a homogeneous random field, (b) the 
assumption of ergodicity of this field was also introduced. 

Weak homogeneity implies that the mean initial imperfection 
function is a constant, whereas the autocorrelation function 
depends only on the difference x 2 - x lf where Xi and x 2 are spatial 

coordinates. This insensitivity to the shift of initial cross- 
section of reference is possible for infinite domains. Therefore, 
possibly, this infinite length assumption that was adopted. For 
solving the problem authors resorted to the classical method of 
stochastic linearization and the additional assumption, that the 
output random field, namely the additional deflection of the 
column, was ergodic too. The main conclusion derived in the paper 
was that each infinite column in the ensemble has the same buckling 
load, which depends on the autocorrelation function of the initial 
imperfection alone, not on a particular realization of any of them. 
Since then, perhaps because of this surprising conclusion, this 
problem attracted the attention of other investigators who tackled 
the problem by various methods. Authors used the method of 
stochastic linearization (Amazigo et al, 1971) , truncated hierarchy 
(Amazigo, 1969, Amazigo et al, 1971), and perturbation (Amazigo, 
1971; Amazigo 1974). In his review Amazigo (1976) comments on the 
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ergodicity assumption: 

“ The method . . .based on the ergodicity hypothesis leads to the 
conclusion that the structure will buckle statically or 
dynamically at the corresponding ( deterministic , I.E.) load... 
with probability 1. This result may appear paradoxical . 
However , to dispel the apparent contradiction we note that no 
matter how the origin of an infinitely long column is defined 
the buckling load for such columns with imperfections 
w(x) = sin(x + <t> ) is independent of <p and hence independent of 
any probabilistic distribution we may assign to 
It appears to us that the source of this paradoxal result 
stems from the fact that the authors assumed the ergodicity not 
only of the input field, but also of the output field. This 
assumption allowed to facilitate the solutions that were derived. 
In order to check the validity of such an assumption Scheurkogel et 
al (1981) undertook an investigation of a model system, which 
allowed to obtain a closed-form solution. Then the same problem 
was solved by invoking the ergodicity assumption. A control 
parameter k was introduced so that one could study the varying 
behavior of the system as the control parameter was changed. It 
turned out that, in general the output of the system was inergodic. 
At some value of the parameter, k = 2, the ergodicity assumption 
yielded a result coinciding with the response obtained exactly. 
This implies that sometimes the error, may not affect the estimate 
of the system's response! In two distinctive ranges of parameter 
k the behavior turned out to be of different nature. ForO < k < 2 
the ergodicity assumption introduced a small error of the order of 
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one percent. Yet, for 2 < k <, 4 the ergodicity assumption 
introduced a small error of the order of one percent. Yet, for 
2 < k £ 4 the ergodicity assumption led to large errors. In 
particular, when k tends to four the ergodicity-based solution is 
finite, whereas the exact solution is unbounded. As is seen, 
extreme caution must be exercised when invoking the ergodicity 
assumption: the differential equation itself, rather than an 
analyst, will decide if the output is ergodic or not! 

Closely related conclusions were arrived at in another 
investigation (Scheurkogel et al, 1985) which studied an applied 
mechanics problem, in which Bolotin (1971) also utilized an 
ergodicity assumption. 

3. MONTE CARLO METHOD 


The present author was introduced to stochastic buckling 
rather incidentally. During the academic year 1977/78, Professor 
Bernard Budiansky was supposed to spend a sabbatical year at the 
Department of Aerospace Engineering of the Technion-Israel 
Institute of Technology. Head of the Structures Group, Professor 
Joseph Singer recommended me to investigate some problems which 
may generate an interest of the sabbatical visitor. Without taking 
an obligation to do so (we all have an academic freedom, don't we?) 
I decided to study some of the works of Professor Budiansky in more 
detail. 

I read several elegant articles . of Budiansky and Hutchinson 
voted to the imperfection sensitivity of structures. Then, when I 
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read the article by Frazer and Budiansky (1969) , their result of it 
struck me: the realizations of the columns are different, yet they 
all share the same, deterministic buckling load! Yes, this load 
did depend on the probabilistic characteristic of the imperfection, 
yet it was a constant! Interestingly and surprisingly, this load 
depended upon the single value of the' spectral density of initial 
imperfections, and was independent, otherwise, on the spectral 
content of the field. 

I have decided that instead of pursuing a new purely analytic 
approach, in addition to what was already undertaken in Harvard, it 
would be nice to perform an experiment. Yet, where should we get 
numerous realizations, I thought, of the real column on nonlinear 
elastic foundations? If not in the real laboratory, then may be in 
the virtual one: on the computer? Thus the idea occurred to study 
the Frazer-Budiansky problem by the Monte Carlo simulation. 

The idea did not seem to be very fancy or even new. Indeed, 
Frazer (1965) performed the Monte Carlo analysis of a column on a 
nonlinear foundation. Yet, most unfortunately, he limited himself 
with a single-term approximation which leads, once this assumption 
is made, to a closed-form solution. Naturally, one does not need 
the Monte Carlo solution if the exact solution is at hand, except 
when one wants to illustrate the validity of the Monte Carlo 
solution in the particular case capable of the exact solution. 
Once the confidence is gained, one resorts if the multi-mode 
solution when the exact solution is unavailable. 

Multi-term Monte Carlo simulation was conducted by Hansen 
(1977) in his probabilistic analysis of randomly imperfect shells. 
However, the analysis performed could be characterized as an 


12 



unbalanced ones fancy analytical analysis was superimposed with a 
naive probabilistic analysis. The assumption was made that all the 
Fourier coefficients used in the series expansion were identically 
distributed. Namely, each Fourier coefficient was taken as a 
normally distributed variable with the same variance. This 
assumption in essence corresponds to the "white-noise" 
autocorrelation function of the initial imperfections. Thus, the 
analysis did not allow the information on general autocovariance 
function. Some other investigators too, although in a dynamic 
context, neglected correlations between the various Fourier 
coefficients, although adopted the nonconstant variances (even a 
new term was coined for this doubtful proposition: "grey noise," 
see Lindberg, 1983) . 

It was realized by this writer that the simulation analysis 
should start from the mean function and the autocovariance 
function, and end up with the variance-covariance matrix of initial 
imperfection's Fourier coefficients. This matrix, in general case, 
must be a fully populated one: not a diagonal one with identical 
(Hansen 1977) or different (Lindberg, 1993) elements. In order to 
study the Frazer-Budiansky model structure the present writer 
developed a general simulation procedure for solving the stochastic 
boundary value problems (Elishakoff 1979) . 

This simulation procedure was applied to the impact buckling 
of a column (Elishakoff 1978) , Hoff's problem of buckling of a 
column in a testing machine (Elishakoff, 1980) and to the Frazer- 
Budiansky problem (Elishakoff, 1979) . A column of finite length 
was studied because of several reasons: (1) an assumption of 
infinite length may simplify the analytical analysis but may 
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complicate the numerical one, (2) the structures utilized and 
analyzed by engineers do not p'osses with the infinite length, (3) 
the nonlinear column on nonlinear elastic foundation does not have 
a behavior of the edge-effect to justify for looking for interior 
solutions as those associated with the infinite structure, and edge 
effect solutions near the edges. For each realization buckling 
load was found numerically by transforming the nonlinear algebraic 
equations to the numerically solved ordinary differential 
equations. The adopted method of Qiria (1951) and Davidenko (1953) 
is similar to the arc-length method of Riks (1979 ) , as was 

communicated to this writer by Professor W. T. Koiter. Reliability 
of the column was calculated. Following conclusions were drawn: 

(1) Monte Carlo solution yields results that are practically 
coincident with the exact solution, when the latter is 
available. This demonstrates that the Monte Carlo solution 
may exhibit a better performance than the various statistical 
tests may predict. 

(2) A single-term Galerkin approximation is not sufficient to 
reliably predict the structural reliability; depending on the 
system's parameters, various but higher degrees of 
approximation must be achieved in order reliability estimates 
to be accurate. 

(3) Design buckling load associated with high reliability may 
significantly deviate from the average buckling load. 

(4) When the length of the column increases the variance of the 
buckling load decreases. 

The latter conclusion is not in disagreement with the result of 
Frazer and Budiansky (1969) that for an infinite column the 
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buckling load is a deterministic quantity. For the realistic 
finite column the buckling load depends on the particular 
realization of the initial imperfection function, which in turn 
depends on the probabilistic measures (mean and covariance 
function) of the initial imperfections. 

In a later analysis Day (1980) showed that in some simple 
cases the ergodicity assumption may even be dispensed with for the 
evaluation of the mean buckling load. Yet the analyses yielding 
mean buckling load alone could hardly be considered practical. 
Each of us may remember various, sometimes entertaining, objections 
to average quantities. Anyway, the knowledge of the average 
buckling load is insufficient for probabilistic design of 
structures undergoing buckling. 

Having a general simulation procedure for initial 
imperfections with given mean and covariance functions pinpointed 
the way of introducing the initial imperfection sensitivity into 
design. It involves three main items: 

(a) Development of accurate deterministic (analytical or 
numerical) tools, for buckling load prediction. 

(b) Compiling extensive experimental information on imperfections, 
boundary conditions, elastic properties, scatter in loads 
etc., in view of deriving mean functions of random fields, and 
their covariance functions, and assessing their distributions. 

(c) Utilization of the Monte Carlo analysis with simulating 
brothers and sisters (but not perfect clowns!) of the 
experimentally measured structures. 

This writer was humbled to read in several publications of Arbocz 
(1991) : 
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"It was not until 1979, when Elishakoff published his 
reliability study ... that a method has been proposed, which 
made it possible to introduce the results of imperfection 
surveys into the analysis 
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SECTIONS THAT REMAIN TO BE WRITTEN 

4 . RELIABILITY OF SHELLS OR HOW THE BEAUTIFUL THEORIES ARE KILLED 
BY UGLY EXPERIMENTAL FACTS 

5. A TALE OF THICKNESS VARIATION 

6. IS PROBABILISTIC THEORY SO GOOD, AFTER ALL? (PROBABILITY IS 
NOT A MAGIC WAND) 

7. CONVEX MODELING OF INITIAL IMPERFECTIONS 

8. OPTIMIZATION AND ANTI -OPTIMIZATION UNDER BUCKLING 
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FUTURE NEEDS 


It appears to this writer that the research should concentrate 

on several directions: 

(1) Accumulation of data for statistical analysis to check the 
nature of the distribution of random initial imperfections, 
elastic moduli, thickness variations, load variations etc. 

(2) Development of techniques of the identification of boundary 
conditions, which may turn out to have a nonuniform nature. 
When limited data is provided, the problem of identification 
may be replaced by establishment of local modifications in 
boundary conditions, during the use of the structure, via 
convex modeling. 

(3) Development of finite element codes in stochastic setting, 
incorporating uncertain imperfections, elastic moduli, 
boundary conditions, thickness variation, and loading 
conditions development of buckling post processors to 
commercially available codes like NASTRAN, ADINA, ALGOR etc. 

(4) Possible establishment of a "NASA-University Buckling 
Institute,” which will establish direct connection between the 
industry and researchers, to deal with relevant problems. 

(5) Possible launching of new NASA Initiative, at least on a 
"small fire," so that relevant buckling research would not 
diminish. 
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